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SUMMARY 


A new computer program is presented for calculating the quasi-steady transonic 
flow past a helicopter rotor blade in hover as well as in forward flight. The pro- 
gram is based on the full potential equations in a blade-attached frame of reference 
and is capable of treating a very general class of rotor blade geometries. Computed 
results show good agreement with available experimental data for both straight- and 
swept-tip blade geometries. 


I. INTRODUCTION 


There is an increasing need to develop advanced computational tools for heli- 
copter rotor aerodynamics research. A major thrust to meet this need is to develop 
a reliable and efficient computer code to predict the transonic flow field over a 
helicopter rotor blade. Several investigators, using small disturbance theory, have 
developed computer codes for calculating such flows. Among them, Caradonna and Isom 
(ref. 1) calculated the flow past a nonlifting, hovering rotor blade. Grant (ref. 2) 
considered the quasi-steady flow over a nonlifting rotor blade in forward flight. 
Caradonna (ref. 3) extended his calculation to the unsteady flow past a nonlifting 
rotor blade in forward flight with simple blade geometry. Finally, Chat tot (ref. 4) 
extended Caradonna T s unsteady, small disturbance code to arbitrary blade geometry. 

Arieli and Tauber (ref. 5) were first to publish a code ROT22 based on full 
potential theory for the quasi-steady flow over a rotor blade. Their approach was to 
modify Jameson and Caughey’s widely used fixed-wing computer code FL022 (ref. 6) to 
the case of a rotating blade. This method solves the full potential equations in a 
sheared parabolic coordinate system so that it is possible to treat the blade geom- 
etry exactly. In early comparisons with ONERA data, it was found that pressure dis- 
tributions on the blade were not well-predicted, particularly in the vicinity of 
swept tips (ref. 7). In the process of verifying a correction to the ROT22 code, it 
was expedient to ^develop a new computer code. This new code was denoted TFAR1 
(Transonic Flow Analysis for Rotors) and was found to be a useful tool in its own 
right. It is presented here as an additional method for analyzing the full poten- 
tial, quasi-steady flow on a rotor blade of arbitrary geometry. 

This new code (TFAR1) contains several new features: (1) a new formulation of 

the problem, (2) the capability of treating cranked blades, (3) the capability of 
predicting flow over the blade at any azimuthal angle, (4) the option to restrict 
calculations to the flow over the tip portion of the blade for computational effi- 
ciency, and (5) the option to obtain better resolution by clustering grid points at 
any selected spanwise station. 

The computer results obtained from this new code were compared with ONERA test 
data for straight- and swept-tip blades. The computed results are presented to show 
that (1) the transonic phenomena that take place on the tip of a rotor blade are 



basically three-dimensional and unsteady; (2) the quasi-steady theory predicts good 
pressure distributions for a rotor blade in flow which is either entirely subsonic 
or subcritical with weak shocks; and (3) the quasi-steady theory is useful for design 
work because it gives good pressure distributions for a straight-tip rotor blade near 
the 90° azimuth where the flow has moderate shock waves. 

This report is Part I of a series of planned publications under the same general 
title, "Transonic Flow Analysis for Rotors." 


II. FLOW EQUATIONS 


The exact flow field around a helicopter rotor blade in forward flight is gen- 
erally acknowledged to be a very complex, unsteady, three-dimensional problem. A 
complete numerical simulation is beyond the state of the art. The flow in the pres- 
ent case is assumed to be inviscid and isentropic. Therefore, a velocity potential, 
<J>, exists for the flow described in a frame of reference which is at rest relative 
to the undisturbed air. In this inertial frame, the complete equation for the veloc- 
ity potential is 

4> tt + [(V4>) 2 ] t + V4> • V [j (V<|>) 2 ] = a 2 V 2 <t> (1) 


where a is the local speed of sound. Bernoulli’s equation, relating a and cj) , is 


♦ t + \ (V*) 2 + 




( 2 ) 


where a^ is the sound speed in the undisturbed air, and y is the specific heat 
ratio which is equal to 1.4 for air. 

If the blade geometry and location are described by s(r;t) = 0, where r is 
the position vector in the inertial frame, then the boundary condition at the blade 
surface is 


S t + V<f> • VS = 0 (3) 

For further analysis it is more convenient to implement this surface boundary condi- 
tion in a moving frame of reference in which the blade location is fixed (fig. 1). 

Let primed variables refer to the inertial frame, F ! , and unprimed variables refer 
to the blade-attached moving frame, F. Suppose that at time, t, the two frames are 
coincident and thaj: F is moving relative to F T with a linear velocity, U, and an 

angular velocity, Q. Then, at time, t, the position vector, r, of a particular fluid 
particle is the same for both frames. If a poinjt, P. is rigidly attached in F ’ , it 
is observed in F to move with velocity V = -(U + Q x r) . Thus, the velocity of 
fluid particle at P in F is | = Vf + V. The rate of change of 4> at P is 
measured by an observer in F as 


<f> t » = + t + V • % 

The potential equation in the moving frame, F, is given by 


(4) 
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<(> + 2V • ^<|> + (V - V) (V - V<f>) + V(J) * V t + 2V<J> • V<J> t 

+ 2 1 ( v • V4>)^4>] • ^<J> + V<(> • V (V<j>) 2 J = a 2 V 2 <f> (5) 

and Bernoulli's equation is 

2 a 2 

♦ t + ^ • V4> + \ ( v *) 2 + (6) 

Let the moving frame, F, be described in a Cartesian coordinates system in which 
x, y, and z represent the chordwise, vertical, and spanwise directions of the blade, 
and whose origin Is at the center of rotation. In the inertial frame, let the 
advance velocity, U, lie in the (x T ,z f ) plane, and let it form the inclination angle, 
a Q , with the negative x'-axis direction, and also let the angular velocity, ft, be 
in the positive y’-axis direction. The velocity, caused by the motion of the 
frame, F, has components of 

V x = ftz + U cos a 0 sin 
V 2 = U sin a Q 

and 

V 3 = -ftx 4* U cos a Q cos 

where ^ is the azimuthal angle of the blade = 180° for forward flight direction 
in the inertial frame) . 


The potential equation in Cartesian coordinates is 

<f>,. + 2q d> +2q<f> +2q<J) 

r tt ^i r xt yt b zt 

= (a 2 - q 2 )d) + (a 2 - q 2 )d> + (a 2 - q 2 )(f> - 2q q <J> - 2q q cj) 

4 1 y xx yy *ZZ 4 l 4 2 Y xy 4 1 4 3 T 


xz 


2q q 4> 

M 2 4 3 T yz 


+ (ft 2 x - 2ftU cos a n cos + (ft 2 z + 2ftU cos a sin ip)(p r 


(7) 


where q 19 q 2 , and q 3 are the velocity components of local fluid particle in the 
moving frame and are specified as 


q = 6 + V 

4 l y x 1 

q = (j) + V 

h 2 Y y 2 

and 

q = (J) + V 

4 3 T Z 3 


This equation is similar to the one presented in reference 5 with the exception of 
the last two terms. 
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Bernoulli’s equation in Cartesian coordinates is 


♦t + Vx + 


VJ 
2 T y 


+ V. 


) + i (cf> 2 + 

z 2 Y x 


<j> + <J) ) 

y z 


Y “ 1 Y - 1 


(8) 


The present study is focused on the three-dimensional effect. For steady calcu- 
lation, all time-dependent terms in equations (7) and (8) are dropped. Note that the 
effect of rotation is still included since it is always present in the transforma- 
tion mapping. This steady case is called quasi-steady in the helicopter literature. 
Several boundary conditions are necessary to complete the boundary value problem. 

In the near field, the flow tangency to the blade is described by the expression 

q • n = 0 

where ifr is the normal unit vector to the blade surface. The wake that is shed from 
the trailing edge is assumed to be a vortex sheet which is a smooth continuation of 
the trailing edge. Across this vortex sheet, the pressure is assumed to be continu- 
ous. The jump in potential determined at the trailing edge of each spanwise profile 
is then assumed to propagate to infinity instantaneously. At the far field, the 
boundary condition can be formulated as a Dirichlet condition where the potential 
vanishes . 


III. MESH SYSTEM 


It is much simpler to include boundary condition in a finite difference calcu- 
lation if the boundary surface is conformal with the coordinate surface. A parabolic 
sheared mesh system that is employed in the present analysis is similar to one that 
is used in previous analyses with fixed wings (ref. 6). The mesh system is generated 
by a series of transformations from the physical space to the computational domain 
(fig. 2). 


First, the shearing transformation 

x = x - x (z)^ 
s 

y = y “ y s (z) > 

z = z J 


(9) 


shears out the blade sweep and dihedral. Here, the point x s (z), y s (z) is the center 
of the circle passing through three points near the leading edge of the profile at 
each spanwise station. Second, the scaling transformation 


x = x/SCAL ^ 
y = y/SCAL 
z = z/SCALZ > 


(10) 


accounts for the scaling between the physical space and the computational domain. 
Third, the square root transformation 
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( 11 ) 


(Xj_ + iY x ) 2 = 2(5 + iy) 


Z i = 2 


maps the entire blade surface to a shallow bmp Y x = S(X 13 Z ± ) near the plane 
Y x = 0. Fourth, the second shearing transformation 


X = X x 

Y = Y x - S(X 1# 



( 12 ) 


reduces the blade surface to a portion of the plane Y = 0. Finally, the stretching 
transformations are introduced to render the computational domain finite. For 
example , 


Y = — b>0,0<a<l 

(1 - Y 2 ) a 

is used to map the planes Y = ±°° to Y = ±1. Similar transformations are used out- 
board of the blade tips in the Z-direction and downstream of the trailing edge in 
the X-direction. At the blade trailing edge, the branch cut in each spanwise plane 
is continued smoothly downstream. This cut will be taken as the location of the 
vortex sheet across which the wake condition is applied. 

The transformations (9)-(12) applied to the steady form of equation (7) yield 
an equation of the form 

A +xx + b *yy + c +zz + 2d +xy + 2E *xz + 2f +yz + Vx + Vy + Vz ■ 0 (13) 

If we introduce the following notation, 

0 = x i~ 

X 

V = X. 

y 

K = -(XgO + y^y) 

T| = x s v ~ y s a 
a = -(S x a + y) 

3 = a - S x y 

C = SCAL/SCALZ 

Y = n - ?s x - cs z 
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U = q x a + q 2 y + q 3 £ 
V = q^ + q 2 B + q g Y 


W 

TT 

X 

e 

K 

e 

6 

L 

M 

N 

P 

and 


^3 

y’v + 

y^x - x;y 

s x x + ¥ 

V - * 

x s° + ^ 

x"a + y"3 

s s 


it - eSCAL 



(fi 2 x - 2fiU cos a Q 

COS 40 

• S( 

(f2 2 z + 2fiU cos a Q 

sin 40 

• S( 

X - s x w 




Q = R - 6 • SCAL 

Then, the coefficients of equation (13) can be written as 

A = U 2 - a 2 (a 2 + y 2 + 5 2 ) 

B = V 2 - a 2 (a 2 + g 2 + Y 2 ) 



D = U 
E = V 
F = W£ 

R 1 = crM + £N + (q 2 - + (q* - a 2 )L + Zq^X^ - Zq^x - Zq^I 1 

R 2 = aM + yN + a 2 [(a 2 + y 2 + ? 2 )S XX + 4 2 S ZZ + 2^?S XZ + Q] - (q 2 - q 2 )A 

- 52s xx - ’S'^ZZ - 2ii V s xz - ' 2 1 i< 1 2 £ + 2 »l1, e + 

and 

R 3 = 5N 

At the blade surface, the tangential flow condition is simply V = 0 in the 
X,Y,Z coordinates. At the far field, the Dirichlet condition <j> = 0 is imposed in 
the present study. For points on the continuation of the singular line outboard of 
the blade tips, where the Jacobian vanishes, the potential equation reduces to the 
Laplace equation 


^XX + ^YY 


= 0 


IV. FINITE DIFFERENCE APPROXIMATION 


The potential equation can be rearranged in the canonical form 

(a 2 - q 2 )<(> + a 2 (V 2 <f> - <f> ) 4- first-order terms = 0 (14) 

s s s s 

where q is the magnitude of the velocity, q, and s is the local flow direction. 
This equation is elliptic for subsonic flow (q < a) and is hyperbolic for supersonic 
flow (q > a) . At subsonic points, central differences are employed to approximate 
all derivatives. At supersonic points, upwind differences are applied to <j> ss of 
the first term of equation (14) whereas central differences are used to approximate 
the rest of the terms of equation (14) . In the X,Y,Z system, 

- 02 +XX + 02 *YY + " 2 *ZZ + 2iS *XY + 2D “+XZ + 2 ™+YZ (15) 

It is essential for rotor flow calculation to apply upwind differences to all the 
derivatives in expression (15) in all three directions according to the sign of U, 

V , and W . 


7 



V. SOLUTION ALGORITHM 


A generalized line relaxation scheme is used to solve the finite difference 
approximations of the flow equations in X,Y,Z system. A typical central differ- 
ence formula for (f)^x 


*xx 


4> 


n+l 

i-ijk 


- (2/a))(j) 


n+l 

ijk 


- 2(1 - 

AX 2 


+ <f> 


n 

i+ijk 


where the superscripts denote the iteration level and m is the relaxation factor. 
Similarly , 




n _ .n n+l i ,n+i 

_ , j H-i , k i+ 1 y j ~ j- ?k i^iyj+i^k 


r XY 4 AX AY 

At supersonic points, for positive U and V, typical upwind differences are 




XX 


0 n+1 n 0 ,n+i n 

hjk hjk h-i ,.jk h-2,jk 

AX 2 


and 


,n+i .n+l ,n+i , ,n+i 

, _ hjk h-i , jk hj-i,k h-i,j-i,k 

V XY AX AY 

The relaxation process can be regarded as an approximation to some artificial time- 
dependent equation if we regard each iteration as representing an advance At, in 
artificial time coordinate (ref. 8). An additional term of the form 

0 At<f> gt = g At[U<j> Xt + Vc() yt + W<D zt ] , g > 0 


has been added to this artificial time-dependent equation to speed up the convergence 
rate of the scheme. The upwind differences are used to approximate the spatial 
derivatives of this term. This implicitly introduces a convective viscosity to the 
equation and the scheme is further stabilized. The resulting linear system for the 

unknown is very large. However, its horizontal lines (j and k constant) are 

J-jk 

decoupled. Each horizontal line can thus be solved by a tridiagonal matrix solver. 


VI. RESULTS AND DISCUSSION 


A typical run consists of 100 relaxation sweeps on each of three different grids 
(a finer grid containing twice as many grid points in each direction of a coarse 
grid). An initial calculation is performed on a coarse grid containing 32 by 6 by 8 
grid points in the X, Y, and Z directions, respectively. The solution is then 
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interpolated onto a medium grid and is used as a starting guess. The process is 
repeated again for the fine grid to get the final solution. A typical run for each 
azimuthal position takes about 40 sec (CPU time) on the NASA Ames Cray 1-S computer. 

Comparisons are made with experimental data from two model rotor blades that 
were tested at ONERA in 1978. The detailed blade geometries, one of which had a 
swept tip, are described in reference 9. Both blades are tapered and have symmetric 
blade sections. The swept-tip blade has a 30° leading edge sweep on the outer 15% 
of the blade (the kink is at r/R = 0.85). Their geometries, figures 3(a) and 3(b), 
are approximated in TFAR1 by the respective geometries, figures 3(c) and 3(d). The 
trailing edge of the approximate blades is sawtoothed. 

The first set of results presented is for the nonlifting straight tip blade at 
a free-stream Mach number of 0.2406 (q^ = 80.4 m/sec) and a tip Mach number of 0.5976 
because of rotation (ftR = 199.7 m/sec). The advance ratio is about y = 0.4. Fig- 
ure 4 compares the calculated and measured surface pressure distributions at three 
different span stations, r/R = 0.85, 0.9, and 0.95 for azimuthal angles from 0° to 
180° at 30° increments. Agreement is good for this case. It is noted that the 
flows are either entirely subsonic or subcritical with a small supersonic zone. In 
other words, when the unsteady effect of the flow is small, the code predicts good 
pressure distributions. 

The second set of results that is presented is for the same straight tip blade 
at a free-stream Mach number of 0.3292 (q^ = 110 m/sec) and at a tip Mach number of 
0.5976 because of rotation (ftR = 199.7 m/sec). The advance ratio is high (y = 0.55). 
Figure 5 compares the calculated and measured surface pressure distributions at the 
same three span stations (r/R = 0.85, 0.9, and 0.95) for azimuthal angles from 0° to 
330° at 30° increments. Overall, agreement is fine for the advancing flow side, and 
is poor for the reverse flow side. The flow fields in the advancing flow side are 
subsonic with moderate or greater zones of embedded supersonic flow. The code pre- 
dicts stronger shock waves in the first quadrant and predicts weaker shock waves in 
the second quadrant when compared with the ONERA data. It should be pointed out that 
the code predicts good pressure distributions near 90° azimuth in spite of the 
unsteady effect that is quite strong there. A comparison between TFAR1 and R0T22 
results (ref. 5) at span station (r/R = 0.9) for the azimuthal angles (tJj - 60, 90, 
and 120) is shown on figures 5(c)-5(e). The differences may be due to the absence 
of terms in the flow equation as we mentioned (eq. 7). 

A similar calculation is performed for the swept blade at a free-stream Mach 
number of 0.3127 (q TO = 105 m/sec), a tip Mach number 0.6288 (QR = 210 m/sec) caused 
by blade rotation and a 0° angle of attack. The advance ratio is y = 0.5. Fig- 
ure 6 shows computed and experimental surface pressure distributions for this case. 
The prediction with TFAR1 in the vicinity of the crank is good. TFAR1 accounted for 
the leading-edge sweep and its effect on the pressure distribution. One of the 
effects of the sweep-back of the blade tip is to delay the shock formation. This 
can be seen from the fact that the code TFAR1 predicts good pressure distributions 
at the 120° azimuth for the 30° swept-tip blade. 


9 



VII. CONCLUSIONS 


A finite difference code, TFAR1 , for predicting quasi-steady transonic flow 
over a helicopter rotor blade was presented. The code solves the second order full- 
potential equation in the moving frame and is suitable for modeling the thickness 
effect of the blade. 

Computed results obtained from this new code have been compared with ONERA data 
for both straight- and swept-tip blades with advance ratios ranging from 0.4 to 0.55. 
Results showed excellent comparisons between quasi-steady computations and experi- 
mental pressure distributions for flow which was entirely subsonic or subsonic with 
a small supersonic zone. Fair correlation between quasi-steady computational and 
experimental pressure distributions were obtained for flow with moderate or greater 
zones of embedded supersonic flow. 

It is concluded that (1) quasi-steady theory can predict good pressure distri- 
butions for flows without any shock or with weak shocks, (2) quasi-steady theory can 
still predict good pressure distributions for a straight-tip blade near 90° azimuth 
for flows having moderate shocks and thus is good for design work, (3) the unsteady 
effect that takes place on the tip of a rotor blade on the advancing side is 
basically caused by the transient shock movement, and (4) an unsteady theory is 
necessary to predict the flow field around a helicopter rotor blade when shocks of 
moderate strength appear. 
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APPENDIX A 


DESCRIPTION OF THE CODE 


The input data deck consists of sequences of pairs of cards. The first card of 
each pair gives the names of the parameters that appear on the data cards that follow. 
All data items are read as floating point numbers in a field of 10 columns, and val- 
ues that represent integer parameters are converted in the program. All the input 
data is immediately printed as output so that it is easy to check the accuracy of the 
input . 

After the flight condition is read in, the blade geometry is defined by giving 
blade section profiles at successive spanwise stations from blade root (near the 
center of revolution) to blade tip. The blade planform and dihedral are determined 
by specifying the chord, the leading edge coordinates, and twist angle at each sec- 
tional profile. After the first airfoil is read in, only the leading edge coordi- 
nates, the chord, and the twist angle are given at the new station if this new 
sectional profile is similar. Otherwise, the input profile should be provided again. 
The blade sections between two given stations are generated by interpolation. The 
program prints the coordinates of the unfolded sectional profiles that are produced 
by the code at the root and at the tip of the blade. They should be inspected to 
see if they are reasonably smooth. 

The program also prints a chart of values of an indicator-IV which shows the 
characteristics of points in the Y = 0 plane. The indicator-IV = 2 indicates a 
point on the blade, IV = 1 indicates a point on the trailing vortex sheet, IV = 0 
indicates a point on the singular line, IV = -1 indicates a point adjacent to the 
edge of the blade on the vortex sheet, and IV = -2 indicates an ordinary point 
beyond the blade or vortex sheet. 

The program next displays the iteration history. The maximum correction to the 
velocity potential and the maximum residual of the difference equation together with 
the i, j, and k location, the relaxation factors, the circulation at the middle 
blade section, and the number of supersonic points are printed at every cycle. 

After a specified maximum number of cycles has been completed, or a convergence 
criterion has been satisfied, the section lift, drag, and moment coefficients are 
printed for each span station and the pressure distribution is printed or displayed 
in a plot as desired. Finally, the characteristics of the blade are printed which 
include the coefficients of lift and form drag that are computed by integrating the 
surface pressure. An estimate of friction drag coefficient may be supplied in the 
input, and this will be included to produce an estimate of the total drag coeffi- 
cient. At the end, additional plots are generated if they are desired. These show 
a view of the blade and the three-dimensional pressure distributions over the upper 
and lower surfaces, respectively, with the root at the bottom of the picture. 
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APPENDIX B 


GLOSSARY OF INPUT PARAMETERS 

TITLE Title of the case being run. (A format) 

Card pair 1: 

FNX The number of mesh intervals in the direction of the chord. 

FNY The number of mesh intervals in the direction normal to the chord and 

span. 

FNZ The number of mesh intervals in the direction of the span. 

Card pair 2 : 

FIT The maximum number of iteration cycles which will be computed. 

COVO The desired accuracy. 

P10 The subsonic relaxation factor for the velocity potential. P10 lies 

between 1 and 2 and should be increased toward 2 with mesh refinement. 

P20 The supersonic relaxation factor for the velocity potential. Recom- 

mended value 1. 

P30 The relaxation factor for the circulation. Recommended value 1. 

BETAO The damping factor which controls the amount of convective term. Recom- 

mended value 0.1. 

FHALF Determines whether the mesh will be refined. FHALF = 0 terminates the 

computation after FIT iterations or after convergence. FIT = 0 
halves the mesh after FIT iterations or convergence on the coarse 
mesh. An additional card pair 2 is required for each mesh refine- 
ment. The value FHALF = 0 appears on the last mesh refinement 
card . 

Card pair 3: 

FSPEED The forward flight speed (m/sec) . 

PSI The azimuthal angle of the blade (deg) . 

ALPHA The angle of attack (deg) . 

TIPWR The tip speed due to the rotation of the blade (m/sec) . 

RADIUS The rotor disk radius (m) . 

AINF The speed of sound of undisturbed air in far field (m/sec) . 
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Card pair 4 

CREF 

XREF 

FBLADE 

FCLUST 

CDO 

Card pair 5 
FNC 
SWEEP 1 
SWEEP 2 
SWEEP 
DIHED1 
DIHED2 
DIHED 

Card pair 6 

ZS 

XL 

YL 

CHORD 

THICK 

TWIST 

FSEC 


The reference chord length. 

The reference chordwise ordinate of point about which the sectional 
airfoil pitching moment coefficient is calculated. 

Controls the tip portion of blade to be calculated. FBLADE = 1 gives 
the whole blade. 

Controls the spanwise mesh point distribution. FCLUST = 0 means uni- 
form grid is used. 

The estimated drag due to skin friction. This can be added to the drag 
calculated by the code to give the total drag. 


The number of span stations from the blade root to the tip. 
The sweep of the singular line at the blade root (deg) . 

The sweep of the singular line at the blade tip (deg). 

The sweep of the singular line in the far field (deg) . 

The flap angle of the singular line at the blade root (deg) . 
The flap angle of the singular line at the blade tip (deg) . 
The flap angle of the singular line in the far field (deg). 


Span location of the section. 

X coordinate of the leading edge. 

Y coordinate of the leading edge. 

The, local chord value by which the profile coordinates are scaled. 

Modified the section thickness. The Y coordinates are multiplied by 
THICK. 

The angle through which a section is rotated to introduce twist about 
the quarter-chord point of the section airfoil. 

Indicates whether or not the geometry for a new profile is supplied. 

FSEC = 0 means the section is obtained by scaling the profile used 
at the previous span section according to the parameters CHORD, THICK, 
AND TWIST. No further cards are read for this span station and the 
next card is the title card for the next span station, if any. 

FSEC = 1 means the coordinates for a new profile are to be read 
from the data cards that follow. 
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Card pair 7: 

YSYM Indicates the type of profile. YSYM = 0 means the data supplied are 

for a cambered profile. Coordinates are given for the upper and the 
lower surfaces, each ordered from nose to tail with the leading edge 
included in both surfaces. YSYM = 1 means the data supplied are 
for a symmetric profile. A table of coordinates is read in for the 
upper surface only. 

FNU The number of upper surface coordinates. 


FNL The number of lower surface coordinates. For YSYM « 1, NL = NU. 

Card pair 8: (Upper surface coordinates) 


X,Y The coordinates of the upper surface. They appear from leading edge to 

trailing edge. 

Card pair 9: (Lower surface coordinates) 

X,Y The coordinates of the lower surface from leading edge to trailing edge. 

The leading edge of the upper surface is the same as the leading edge 
of the lower surface. The trailing edge points are different if the 
profile has an open tail. 

Card pairs 10, 11, . . .: 

These card pairs are like card pairs 6, 7, 8, and 9. The number of such 
card pairs depends on the number of span stations, FNC. 
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APPENDIX C 


LISTING OF TFAR1 PROGRAM 


♦COMDECK BLANK 
COMMON/ / 


♦COMDECK A 

COMMON/A/ 


♦COMDECK FLO 


G(129,18,33),S0(129,33),E0(33), 

IV (129 ,33) , ITE1(33),ITE2(33),SMACH(33) , 

AO (129) ,A1 (129) ,A2(129) , A3 (129) , 

BO ( 1 8 ) ,B1(18) ,82(18) ,B3(18), 

CO (33) , Cl (33) ,C2(33) ,C3(33) , 

XC(33),XZ(33),XZZ(33),YC(33),YZ(33),YZZ(33), 

NX,NY,NZ,KTE1,KTE2,1SYM,KSYM,SCAL,SCALZ, 

ALPH A, FMACH, HINGE, OMEGA, PSI,CA,SA,CAC, CAS, 
TMACH , TILT , RAD , PI , NIT , PSIS , 
AAO,K2,K3,LX,MX,MY,MZ,DX, DXX , WATX , WATY 

SX(129,33),SZ(129,33), 

SXX(129,33) ,SXZ(1 29,33 ) ,SZZ( 129, 33) 


COMMON/FLO/ PI , P2 , P3 , BETA , FR , IR , JR , KR , GD , IG , JG , KG , NS 


♦DECK TFAR1 

PROGRAM TFAR1 ( INPUT , OUTPUT , TAPE5=INPUT , TAPE6=OUTPUT , 

, TAPE14,TAPE16,TAPF.20) 

C THREE DIMENSIONAL ROTOR BLADE ANALYSIS IN TRANSONIC FLOW 

C USING SHEARED PARABOLIC COORDINATES 

C G IS THE REDUCED VELOCITY POTENTIAL IN THE MOVING FRAME 

C A ROTATIONAL FLOW VERSION OF JAMESON'S FL017 FOR A ROTOR BLADE 

C PROGRAMMED BY I-CHUNG CHANG, APRIL 1983 

♦CALL BLANK 

♦CALL A 

♦CALL FLO 

DIMENSION XS(200,33),YS(200,33), 

1 ZS(33) ,XLE(33) ,YLE(33) , SLOPT ( 33 ) , TRAIL ( 3 3 ) , NP ( 33 ) , 

2 El (33) ,E2(33) ,E3(33) ,E4(33) ,E5(33) , 

3 XP(200) ,YP(200) ,D1 (200) ,D2(200) ,D3(200) , 

4 X(129) ,Y(129) , SU( 129) ,SV ( 129) , SW ( 1 29) ,SM(129) , 

5 CP (129) , CHORD (33) ,SCL(33) ,SCD(33) ,SCM(33) , 

6 FIT ( 3 ) , COV 0 (3) ,P10(3),P20(3) ,P30(3) ,BETA0(3) , 

7 FHALF(3) ,FITMIN(3) , K PLOTS (26),TYTLE(10) 

ND = 200 

NE = 129 

IREAD a 5 

IWRIT = 6 

PI = 3.141592b5358979 

RAD a 57.2957795130823 

DRAD a 1 . /RAD 
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1 WRITE ( IWRIT , 600 ) 

WRITE (IWRIT, 2) 

2 FORMAT ( 1 4HOPROGR A-M TFAR1 , 70X , 3 1 H 1-CHUNG CHANG , NASA-AMES CENTER/ 

1 42H0THREE DIMENSIONAL ROTOR BLADE ANALYSIS IN, 

2 51H TRANSONIC FLOW USING SHEARED PARABOLIC COORDINATES) 
READ(IREAD,530) TYTLE 

WRITE(IWRIT,b30) TYTLE 
READ(IREAD,500) 

REA0(IREAD,510) FNX , FN Y , FNZ 
NX = FNX 

NY = FN Y 

NZ = FNZ 

IF (NX.LT.l) Gu TO 302 
I PLOT = -1 

KPLUT = 1 

READ ( I RE AD, 500) 

NM =0 

11 NM — N M +1 

READ(IREAD,510) FIT(NM) ,COVO(NM) , PI 0 ( NM ) , P20 ( NM ) , 

. P30(NM) ,BETA0(NM) ,F HALF(NM) 

IF (FHALF(NM) .NE.O. .AND.NM.LT. 3) GU TO 11 
FH AbF ( 3 ) = 0. 

READ(iREAU,500) 

READ ( I READ , 510) FSPEED , PS I , ALPHA , TI PWR , RADIUS , A INF 
READ(IKEAD,500) 

READ (IREAD, 510) CHEF , XREF , FBLADE , FCLUST , CDO 

TILT = ALPHA 

ICLUST = FCLUST 

ALPHA = TILT ♦DRAD 

PSID = PSI 

PSI = PSID *DR AD 

DRADIUS s 1. /RADIUS 

FMACH = FSPEED/ A INF 

TMACH = T I P * R / A I N F 

POMEGA s riPwR/(CKEF*KADIUS) 

UTKEF = CREF/AINF 

OMEGA = POMEGA+UTREF 

CA s FMACH*C0S( ALPHA) 

SA = FMACH*SIN( ALPHA) 

PSIM a PSI-.5*PI 

CAC = C A ♦COS ( PSIM ) 

CAS = C A*S1 N ( PSIM ) 

CALL GEOM ( ND , NC , NP , ZS , XS , YS , XLE , YLE , SLOPT , TRAIL , XP , YP , 

1 SWEEP1 ,SWEEP2, SWEEP, D1HED1 ,DIHED2 ,DIHED, 

2 D1 ,D2,XTEO,CHORDO, ZTIP, IS YMO, HINGE, FBLADE) 

IS YM = IS YMO 

91 CALL COORD ( NX , N Y , NZ , XTEO , ZTIP , XMAX , ZM AX , ICLUST , 

1 SY,SCAL,SCALZ, AX, AY ,AZ, 

2 A0,A1,A2,A3,B0,B1,B2,B3,C0,C1,C2,C3) 

CALL SINGL ( NC , N Z , KTE1 , K f E2 , CHGRDO , 

1 SWEEP 1 ,SWEEP2, SWEEP, DIHED1 , DIHED2 , DIHED , 

2 ZS,XLE, YLE,XC,XZ,XZZ, YC,YZ, YZZ, 

3 CO, Cl ,C2,C3,E1 ,E2,E3,E4,E5,IND) 

CALL SURF (ND,NC,ZS, SLOPT, TRAIL, XS,YS,NP, 

1 XP,YP,D1,D2,D3,X,Y,IND) 
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IF (IND.EQ.O) GO TO 291 
NM = 1 

NIT = 0 

CALL ESTIM 

101 WRITE (IWRIT , 600 ) 

MIT = FIT(NM) 

KIT = MIT 

IF (NM.GT.l.AND.FHALF(NM) .EQ.O.) KIT = 10 
JIT = NIT 

KRES = (MIT -NIT -2)/500 +2 

JRES = 0 

NRES * 0 

COV * COVO(NM) 

COU = 10000000. 

KY a NY +1 

K 1 = 2 

K2 = NZ 

103 LZ = NZ/2 +1 

WRITE (IWRIT, 104) 

104 FORMAT (49H01N DICAT ION OF LOCATION OF BLADE AND VORTEX SHEET, 

1 27H IN COORDINATE PLANE Y a 0./ 

2 27H0 ( ( I V ( I , K ) , K=K 1 ,K2) ,1 = 2, NX) ) 

DO 106 1 = 2, NX 

106 WRITE (IWRIT, 650) ( I V ( I , K ) , K=K1 , K 2 ) 

WRITE (IWRIT, 600) 

WRITE (IwRIT,112) 

112 FORMAT (49HOCHORDW1SE CELL DISTRIBUTION IN SQUARE ROOT PLANE, 

1 54H AND MAPPED SURFACE COORDINATES AT CENTER LINE AND TIP/ 

2 15H0 AO , 1 5H A1 , 

3 15H A2 , 1 5H A3 , 

3 15H ROOT PROFILE, 15H TIP PROFILE ) 

DO 114 1=2, NX 

114 WRITE (IWRIT, 610) AO ( 1 ) , A 1 ( I ) , A2 ( 1 ) , A3 ( I ) , SO ( I , KTE1 ) , SO ( I , KTE2 ) 
WRITE (IWRIT, 116) 

116 FORMAT ( 1 5H0 TE LOCATION ,15H POWER LAW ) 

WRITE (IWRIT, 610) XMAX,AX 
WRITE (IWRIT, 600) 

WRITE (IWRIT, 118) 

118 FORMAT (46H0N0RM A L CELL DISTRIBUTION IN SQUARE ROOT PLANE/ 

1 1 5H0 BO , 1 5H B1 , 

2 15H B2 , 1 5H B3 ) 

DO 120 J=2 , K Y 

120 WRITE (IWRIT, 610) B0(J) ,B1(J),B2(J) ,B3(J) 

WRITE (IWRIT, 122) 

122 FORMAT ( 1 5H0 SCALE FACTOR, 15H POWER LAW ) 

WRITE (IWRIT, 610) S Y , A Y 
WRITE (IWRIT, 600) 

WRITE (IWRIT, 124) 

124 FORMAT (27H0SPANWISE CELL DISTRIBUTION/ 

1 15H0 CO , 1 5H Cl ,15H C2 , 

2 15H C3 ) 

DO 126 K=K1,K2 

126 WRITE (IWRIT, 610) CO ( K ) , Cl ( K ) , C2 ( K ) , C3 ( K ) 

WRITE (IWRIT, 128) 

128 FORMAT ( 1 5H0 TIP LOCATION, 15H POWER LAW ) 
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WRITS (IWRIT,610) ZMAX,AZ 
WRITS (IWRIT’,600) 

WRITE (IWRIT,125) 

125 FORMAT! 14H0SINGULAR LINE/ 

1 1 5H0 X SING , 1 5H Y SING , 

2 15H XZ ,15H YZ ,15H XZZ 

3 15H YZZ ) 

DO 127 K=K1,K2 

127 WRITE (IWRIT,610) XC t K ) , YC ( K ) , XZ ( K ) , YZ ( K ) , XZZ ( K ) , YZZ ( K ) 

WRITE (IWRIT,600) 

WRITE ( I W R IT , 1 3 2 ) 

132 FORMAT ( 35H0IT ERAT I VE SOLUTION — - STEADY MODE) 

WRITE (IwRlT,134) 

134 FORMAT ( 15H0 NX ,15H NY ,15H NZ 

WRITE (IWRIT,640) NX,NY,NZ 
CALL SECOND ( T ) 

WRITE (I*RIT,700) T 
WRITE (I WRIT, 13b) 

136 FORMAT ( 1 5H0 TMACH NO ,15H FMACH NO ,15H TITL ANG 

1 15H AZIMUTHAL ANG ) 

WRITE (IwRIT,610) TMACH , FMACH , TILT , PSID 
WRITE ( IWKIT, 138) 

138 FORMAT ( 1 OHO ITERATION # 1 5H CORRECTION ,4H I ,4H J ,4H K , 

1 15H RESIDUAL ,4H I ,4H J ,4H K , 

2 I OH CIKCULATN,lOH REL FCT 1,10H REL FCT 2,10H REL FCT 3 

3 10H bETA ,10H SONIC PTS) 

141 NIT = NIT +1 

JIT = JIT +1 

PI = PIO(NM) 

P2 = P20 ( NM ) 

P3 = P30(NM) 

BETA = BKTAO(NM) 

CALL RELAX 


WRITE ( IWKIT , 660 ) N IT r GD , IG , JG , KG , FR , IR , JR , KR , EO ( LZ ) , 
1 P1,P2,P3,BETA,NS 

IF (.NIT.LT.MIT.AND.GD.GT.COV.AND.GD.LT. 10. ) GO TO 141 
IF (FHALF(NM) . EQ. 0 . ) GO TO 176 



) 
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WRIT£CIWRIT,60Q) 

WRITE(IWRIT, 184) PSID 

184 FORMA T(lrtO,*AZlKUTHAL ANGLE * * , FI 5 . 5 ) 

171 K s K +1 

IF (K.EO.MZ) GO TO 191 

IF (K.LT.KTE1 .OR.K.GT.KTE2) GO TO 171 

11 = ITEl(K) 

12 = IT E 2 C K ) 

ZSEC a COCK) + HINGE 

VROT = OK EGA *Z SEC 

VTAN = VROT + FMACH*SIN(PSI)*COS(ALPHA) 

SMACH(K) = VTAN 

CALL VELO(K,SU,SV , SW , SM , CP , X , Y ) 

175 CHORD(K) r X(I1) -XCLX) 

CALL FORCFCI1,I2,X,Y,CP,TILT,CHORDCK) ,XCCK) ,SCLCK) ,SCDCK) ,SCMCK)) 
CALL PSURE(1PL0T,K, X, Y, CP , I 1 , 1 2 , SCL C K ) , SCD(K) , SCM(K) ) 

WRITE (IWRIT,600) 

WRITE ( IWRIT , 182 ) 

182 FORMAT (24H0S EC T ION CHARACTERISTICS/ 

1 15H0 SPAN STATION, 15H CL ,15H CO , 

2 15H CM ) 

ZPHYS= COCK) + HINGE 

185 WRITE CIWRIT,610)ZPHYS,SCLCK) ,SCDCK) ,SCMCK) 

IF CKPLOT.GE.O) CALL CPLOT Cll ,12,SMACHCKj ,X, Y,SU,SV,SW,SM,CP) 

GO TO 171 

191 CALL TOTFORCKTEl,KTE2, CHORD, SCL, SCO, SCM, CO, XC, 

1 CL,CD1 ,CMP,CMR,CMY) 

CD = COO +CD1 

VLD1 = 0. 

IF CABSCC0D.GT.1.E-6) VLD1 = CL/CD1 
VLD = 0. 

IF CABSCCD).GT.l.E-6) VLD = CL/CD 
WRITE CIWRIT,600) 

WRITE CIWR1T,192) 

192 FORMAT C22H0BLADE CHAR ACTERISTICS/ 

1 1 5H0 CL ,15H CD FORM ,15H CD FRICTION , 

2 15H CD ,15H L/D FORM ,15H L/O ) 

WRITE CIWRIT,610) CL , CD 1 , CDO , CD , VLD1 , VLD 

WRITE (IWRIT 196) 

196 FORMAT C 1 5H0 * CM PITCH ,15H CM ROLL ,15H CM YAW ) 

WRITE Cl WRIT, 610) CMP,CMR,CMY 
210 CALL THREED CIP LOT, SU,SV,SW,SM, CP, X,Y,TYTLE, CHORD, 

. CL,CD,CHORD0,SCL,SCD,SCM) 

CALL ROTORB C I PLOT ,SU,SV,S«,SM,CP,X,Y) 

GO TO 301 

291 WRITE Cl WRIT, 600) 

WRITE C IWRIT , 292 ) 

292 FORMAT C24H0BAD DATA, SPLINE FAILURE) 

301 IFCIPLOT.EO.O) CALL DONEPL 

302 STOP 

500 FORMAT C 1 X ) 

510 FORMAT C8F10.6) 

511 FORMAT C 261 3 ) 

530 FORMAT C 1 0A8 ) 

600 FORMAT C 1 H 1 ) 
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610 FORM AT ( F 12.4,7 K 15. 4) 

620 FORMAT ( 8E1 5 . 5 ) 

630 FORMAT(1HO, 10A8) 

640 FORMAT (18, 7115) 

650 FORMAT (IX, 3214) 

660 FORMA T (110, E15. 5, 314, £15.5,314, 5 F 10. 5, 110) 

661 F0RMAT(I10,E15. 5, 314, £15.5,314, 110) 

670 FORMAT (2E15.4,2F15.4) 

700 FORMA T( 15H0C0MPUTING TI ME , FI 0 . 3 , 1 OH SECONDS) 

900 FORMAT ( 1 X , F 1 2 . 8 , 5X , FI 2 . 8 , 5X , FI 2 . 8 , 5X , F’l 2 • 8 , 5X , F 1 2 • 8 , 5X , 
. F12.R,5X,F12.8) 

END 


♦DECK GEUM 

SUBROUTINE GEOM ( ND , NC , NR , ZS , XS , X S , XLE , YLE , SLOPT , TR AI L , XP , YP , 

1 SWEEP1 ,SWEEP2, SWEEP, DIHED1 , DI HED2 , DIHED , 

2 D1 ,D2,XTEO,CHOKDO, ZTIP,ISYMO, HINGE, FBLADE) 

C GEOMETRIC DEFINITION OF ROTOR BLADE 

DIMENSION XS(ND,1),YS(ND,1),ZS(1),XLE(1),YLE(1),D1C1),D2(1), 

1 SLOPT ( 1 ) , TR Alb ( 1 ) , XP ( 1 ) , YP ( 1 ) , NP ( 1 ) 

IREAD = 5 

IWR1T = 6 

RAD = 57.2957795130823 

READ(IREAD,500) 

READ( IREAD, 510) FNC , SWEEP 1 , SWEEP2 , SWEEP , DIHED1 , DIHED2 , DIHED 
IF (FNC.LT.3.) RETURN 
NC = FNC 

WRITE (IWRIT,2) 

2 FORMAT ( 1 5H0 SWEEP(l) ,15H SWEEP(2) ,15H FINAL SWEEP , 

1 1 5H DIHED ( 1 ) , 1 5H DIHED(2) ,15H FINAL DIHED ) 

WRITE (IWRIT,610) SWEEP1 , SWEEP2 , SWEEP , DIHED1 , DIHED2 , DIHED 


SWEEP1 


SwEEPl/RAD 

SWEEP2 

= 

SWEEP2/RAD 

SWEEP 

= 

SWEEP/RAD 

D1HED1 

s 

DIHED1/RAD 

DIHED2 

s 

DIHED2/RAD 

DIHED 

r 

DIHED/RAD 

ISYMO 

= 

1 

XTEO 

ZL 

0 . 

CHORDO 

2 

0 . 

K 

s 

1 


11 READ(IREAO,500) 

READ (IREAD, 510) ZS ( K ) , XL , YL , CHORD , TH ICK , TWIST , FSEC 

AL = TWIST 

ALPHA = AL/KAD 

IF (FSEC.EQ.O.) GO TO 31 

READ (IREAD, 500) 

READ (IREAD, 510) YSYM,FNU,FNL 
NU = FNU 

NL = FNL 
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N = NU +NL 

READ ( I READ , 500 ) 

DO 12 I=NL,N 
12 READ C IREAD , 51 0 ) XPCI),YPCI) 

L = NL +1 

IF CYSYM.GT.O.) GO TO 15 
READ ( I READ , 500 ) 

DO 8 1 = 1 ,NL 

READ ( IREAD, 510) V A L , D U M 
J = L -I 

XPCJ) = VAL 

8 YP C J ) = DUM 

GO TO 21 

15 J = L 

DO 16 I=NL , N 

J = J -1 

XPCJ) = XP(I) 

16 YP(J) = - YP C I ) 

21 WRITE (IWRIT,600) 

WRITE C IWRIT , 22 ) ZS(K) 

22 FORMAT (16riOPROFILE AT Z = ,F10.5/ 

1 1 5H0 TE ANGLE ,15H TE SLOPE , 15H X SING , 

2 1 5ri Y SING ) 

CALL SINGPTCXP,YP,NL,N,XSING,YSING,TRL,SLT) 

WRITE ClWRIT,blO) TRL,SLT,X5ING,VSING 
WRITE (IWRIT,24) 

24 FORMAT ( 1 5H0 X ,15H Y ) 

DO 26 1=1 , N 

26 WRITE C I WRIT , 6 1 0 ) XPCI),YPCI) 

31 SCALE = CHORD/ C XP Cl) -XPCNL)) 

DO 33 1=1 ,N 

D1C1) = 'XL + SCALE*CXPQ)-XPCNL) ) 

33 D2CI) = YL ♦ SCALE*CYPCI)-YPCNL) )*THICK 

CALL SINGPTCD1 , D2 , NL , N , XSiNG , YSING , TRL , SLT ) 

XLECK) = XSING 
YLECK) = YSING 
CA = COSCALPHA) 

SA = SINCALPHA) 

DO 32 1=1, N 

XSCI,K) = CD1CI) -XSING)*CA +CD2CD -YSING)*SA 

32 YS C I , K ) = CD2CI) -YSING)*CA -CD1CI) -XSING)*SA 

SLOPTCK) = SLT -TANCALPHA) 

TRAILCK) = TRL/RAD 
NPCK) = N 

XTEO = AMAX1CXTE0,XSC1,K)) 

CHORDO = AKAX1 CCHORDO, CHORD) 

WRITE C IWRIT , 52 ) ZSCK) 

52 FORMAT C 27 HOSECTION DEFINITION AT Z = ,F10.5/ 

1 1 5H0 XLE , 1 5H YLE ,15H CHORD 

2 1 5HTHICKNESS RATIO, 15H TWIST ANGLE ) 

WRITE Cl WRIT, 610) XL , YL , CHORD , TH ICK , AL 

K = K +1 

IF CK.LE.NC) GO TO 11 

65 ZO = C1.-.5*FBLADE)*CZSCNC)-ZSC1)) +ZSC1) 

KK =0 
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ZTIP = ZS(NC) -ZO 

DU 63 K=1,NC 

ZS ( K ) = ZS ( K ) -ZO 

IF(ABSCZSCK) ) .GT.ZTIP) GO TO 63 

KK = KK +1 

N = NP(K) 

DO 64 1=1, N 
XS(I,KK) = XS ( I , K ) 

YS(I,KK) = YS ( I , K ) 

64 CONTINUE 

ZS(KK) s ZSCK) 

XLE(KK) = XLE(K) 

YLE(KK) = YLKCK) 

SLUPTCKK) = SL.OPTCK) 

TRAIL C KK ) = TRAIL CK) 

NP(KK) = NPCK) 

63 CONTINUE 

NC = KK 

HINGE = ZO 

RETURN 
RETURN 

500 FORMAT(lX) 

510 FORMATC3F10.6) 

600 FORMAT ( 1 H 1 ) 

610 FORMAT(F12.4,7F15.4) 

END 


♦DECK COORD 

SUBROUTINE COORD ( N X , N Y , NZ , XTEO , ZTI P , XMAX , ZM AX , ICLUST , 

1 SY,SCAL,SCALZ,AX,AY,AZ, 

2 A0,A1 ,A2,A3,B0,B1 ,B2,B3,C0,C1 , C2 , C3 ) 

C SETS UP STRETCHED PARABOLIC AND SPANWISE COORDINATES 

DIMENSION A0Cl),AlCt),A2(l),A3(l),B0Cl),Bl(l) , B2 ( 1 ) , B3 ( 1 ) , 
1 C0C1) ,C1 (1) ,C2(1) ,C3C1) 

PI = 3.14159265358979 

DX = 2. /NX 

DY = 1./N1 

DZ = 2./NZ 

DDX = l./DX 

DDXX = DDX* DDX 

DDY = l./DY 

DDZ = l./DZ 

KY = NY +1 

AX = .5 

AY = .5 

AZ = .5 

XMAX = .625 

ZM AX = .625 

SY = .5 

SCAL = XTEO/ ( .50001*XMAX*XMAX) 
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13 


12 


22 


30 


40 


SCALZ 

= ZTIP/(1.000001*ZMAX) 

Ml 

= SCAL/SCALZ 

U2 

= 1 

V2 

= (DX*DDY)**2 

W2 

a (DX*wl*DDZ)*»2 

DO 12 

1 = 2 # NX 

DD 

= Cl -1 )*DX -1. 

8 

= 1. 

IF (A8SCDD) .GT.XMAX) GO TO 13 

DO 

= DD 

D1 

= 1. 

D2 

= 0. 

GO TO 

8 

IF C DD 

.LT.O.) B = -1. 

A 

= 1 ,-(DD-B*XMAX)**2 

C 

= A**AX 

D 

= (AX + AX -l.)*Cl. -A) 

DO 

= b*XMAX+(DD-B*XNAX)/C 

D 1 

= A*C/C1.+D) 

D2 

= -2.*AX*CDD-B*XMAX)*(3.+! 

AO ( I ) 

= DO 

A 1 ( I ) 

= .5*Di*DDX 

A 2 ( I ) 

= D1*D1*U2 

A3 (I) 

= ,5*DX*02 

DO 22 

J=2,KY 

DD 

= (J-2) *DY 

A 

= 1. -DD*DD 

C 

= A** AY 

D 

= (AY +AY -l.)*(l. -A) 

D1 

= A*C/((l. +D)*SY) 

BOCJ) 

= SY *DD/C 

B1(J) 

= ,5*D1*DDY 

82 (J) 

= D1*D1*V2 

B3(J) 

a -AY*DD*DY*(3. +D)/((1. 

IF(ICLUST.EQ.O) GO TO 30 

AH 

= .049 

BH 

= AH/7. 

CH 

= 8.*P1 

DH 

= Pi/7. 

EH 

= 8.*DH 

DO 32 

K=2 r NZ 

DD 

- (K -1 ) *DZ -1. 

B 

= 1. 

IF (ABS(DD) .GT.ZtoAX) GO TO 33 

IF(ICLUST.NE.O) GO TO 40 

DO 

= DD 

D1 

= 1. 

D2 

= o. 

GO TO 

34 

DD 

= . 8* ( DD + ZMAX) 

IF ( DD 

.GT..125) GO TO 45 

A 

= CH*DD 

B 

= BH*SIN ( A ) 

DO 

= DD -B 

D1 

= 1 . / ( 1 . - CH*BH*COS ( A ) ) 


)/((!. +D)*A) 


+D)*A) 
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D2 = 

GO TO 46 

45 A = 

a = 

do = 

D 1 = 

D2 = 

46 DO = 

GO TO 34 

33 IF (DO.LT.O 

A = 

C = 

D = 

DO s 

D1 = 

D2 = 

34 CO ( K) = 

Cl (K) = 

C2 ( K 1 = 

32 C3(N) = 

RETURN 
END 


(b.*DD - l.)*DH 
AH + SIN (A) 

DD + tS 

1./C1.+ AH*EM»COS(A) ) 

D1*EH*EH*B 
1 .25*D0-ZMAX 

) B = -1. 

l.-(DD-B*ZMAX)»»2 

A^AZ 

(AZ +AZ -l.)*Cl. -A) 

B+ZMAX+CDD-B^ZMAXl/C 

A»C/C1.+D) 

-2.*AZ^(DD-b*ZMAX)»C3.+D)/( (l.+D)»A) 

SCAL,Z*DO 

,54D1*W1*DDZ 

Dl*Dl*w2 

,54DZ*D2 


♦DECK S1NGL 

SUBROUTINE SINGL ( NC , NZ , KTE1 , KTE2 , CHORDO , 

1 SWEEPl,SwEEP2,SREEP,DIHEDl,DIHED2 f DIHED / 

2 ZS,XLE, YLE,XC,XZ,XZZ, YC,YZ,YZZ, 

3 C0,Cl,C2 f C3 f El,F2,E3,E4,E5,IND) 

C GENERATES SINGULAR LINE FOR SQUARE ROOT TRANSFORMATION 

DIMENSION ZS( 1) ,XLE(1 ) , YLE( 1 ) ,XC( 1 ) ,XZ( 1 ) ,XZZ(1 ) , 

1 YCC1) , YZC1) , YZZC1) ,C0(1) ,C1 (1) ,C2(1) ,C3C1) , 

2 Elll) ,E2(1) ,E3C1) ,E4(1) , E5(l) 

DO 2 K=1,NC 

E4 ( K ) = 0. 

2 E5(K) = 0. 

K2 = NZ 


11 DO 12 K=2,K2 

IF (CO(K).LT.ZS(l)) KTE1 = K +1 
IF (COCKJ .LE.ZS(NCJ) KTE2 = K 

12 CONTINUE 

B = CHURDO 

51 * TAN (SrtEEPl ) 

52 = TAN(SwEEP2) 

T 1 = TAiMCQlHFDl ) 

T2 = TAN (UIHEU2) 

CALL SPL1F (l,i'‘C,ZS,XLE,El,E2,E3,l,Sl,l,S2,0,0.,IND) 
CALL INTPL CKTE1,KTF.2,C0,XC, 1 , NC , ZS , XLE , El , E2 , E3 , 0 ) 
CAUL INTPL (KTE1,KT£2,C0,XZ,1,NC,ZS,E1,E2,E3,E4,0) 
CALL INTPL (KTE1 , KTE2 , CO , XZZ , 1 , MC , ZS , E2 , E3 , E4 , E5 , 0 ) 
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n n 


CALL SPLIF (1,NC,ZS,YLE,E1,E2,E3,1,T1,1,T2,0,0.,INU) 
CALL IWTPL (KTE1,KTE2,C0,YC,1,NC,ZS,YLE,E1,E2,E3,0) 
CALL INTPL (KTE1,KTE2,C0,YZ,1,NC,2S,E1,E2,E3,E4,0) 
CALL INTPL (KTE1,KTE2,C0,YZZ,1,NC,ZS,E2,E3,E4,E5,0) 

5 = 6+ TAN (SWEEP) 

51 = B + Sl 

52 s B*S2 

T = B*TAN(DIHKD) 

T 1 = B*T 1 

T2 = B*T2 

N = KTE 1 -1 

DO 22 K*2,N 

ZZ = (COCK) -C0(KTE1))/B 

A = EXP ( ZZ ) 

XC(K) = XC(KTEl) + S + ZZ -(SI -S)+(l. -A) 

YC ( K ) = YC(KTEl) + T + ZZ -(Tl -T)»(l. -A) 

XZ(K) = (S +(S1 -S)*A)/B 

YZ(K) s (T + ( T 1 -T)+A)/B 

XZZ(K) = (SI -S)*A/(B*b) 

22 YZZ(K) = (Tl -T)*A/(B*B) 

31 N a KTE2 +1 

DO 32 K=N,K2 

ZZ a ( CO ( K ) -C0(KTE2))/B 

A a EXP ( -ZZ ) 

XC(K ) = XC ( KTE2 ) +S + ZZ + (S2 -S)*C1. -A) 

YC (K ) = Y C ( K T E 2 ) +T + ZZ +(12 -T)*(l. -A) 

XZ(K) = (S +CS2 -S)*A)/B 

YZ(K) a (T + ( T2 -T)»A)/b 

XZZ(K) = -(S2 -S)*A/(B^B) 

32 YZZ(K) a -(T2 -T)*A/(B^B) 

RETURN 

END 


♦DECK SURF 

SUBROUTINE SURF ( ND , NC , ZS , SLOP! , TRAIL , XS , YS , N P , 

1 XP,YP,D1,D2,D3,X, Y,IND) 

INTERPOLATES MAPPED WING SURFACE AT MESH POINTS 
INTERPOLATION IS LINEAR IN PHYSICAL PLANE 
♦CALL BLANK 
♦ CALL' A 

DIMENSION XS(ND,1),YS(NU,1),ZS(1), SLQPT (1),TRAIL(1),X(1),Y(1), 
1 XP(1) ,YP(1 ) ,D1(1) ,D2(1) ,D3(1) ,NP(1) 

PI = 3.14159265358979 

DX = 2. /NX 

LX = NX/2 +1 

MX = NX +1 

MZ = NZ +1 

IVO a 1 

1V1 = -1 

DO 2 K=1 , MZ 
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1TE1(K) 

s 

MX 


ITE2(K) 

= 

MX 


DO 2 1*1, 

MX 



I V ( I , K ) 


-2 

2 

SO ( I , K ) 

= 

0. 


K 

= 

KTE1 


K2 

s 

1 

21 

K2 

= 

K2 +1 


K 1 

= 

K2 -1 


R2 

= 

1. 


IK (ZS(K2) 

"CO ( K ) ) 21,25,23 

23 

R2 

= 

(COCK) -ZS(K1) )/(ZS(K2) -ZS(K1)) 

25 

R 1 

= 

1 . -K2 


C 

= 

R 1 *XS ( 1 , K 1 ) +R2*XS(1,K2) 


CC 

s 

SQKTCCC +C)/SCAL) 


DO 32 1=2 

# NX 


IF U0(l) 

• 

LT.-CC) 11 = I +1 


IF ( AO ( I ) 

• 

LF.CC) 12 = I 

32 

CONTINUE 




1TE1 (K) 

= 

11 


ITE2(K) 

= 

12 


CC 

= 

AO (1 2 ) /CC 


KK 

r 

K 1 


P 

s 

R 1 

41 

N 

= 

N P ( KK ) 


Q 


SORT (XS(1,KK)/C) /CC 


DO 42 1=2 

,NX 

42 

X(i) 

= 

<J*A0(1) 


ANGL 

s 

PI +PI 


U 

= 

1. 


V 

= 

0. 


DO 44 1=1 

> 



R 

= 

SURT(XS(1,KK)**2 + YS ( I , KK ) 4*2 ) 


IF (R.EG. 

0.) 

GO TO 45 


ANGL 

s 

ANGL +ATAN2( (U*YS(1,KK) -V*XS(I, 

1 


( (J*XS ( I , KK ) +V*YS( I , 


U 

= 

XS(I,KK) 


V 

s 

YSCI,KK) 


R 

35 

SORT ( ( R +R)/SCAL) 


XPU) 

= 

R*COS( ,5*ANGL) 


mi) 


R*SIN( ,5»ANGL) 


GO TO 44 



45 

ANGL 

= 

PI 


U 

35 

-1. 


V 

= 

0. 


XPU) 

= 

0. 


YP(l) 

S 

0. 

44 

CONTINUE 




ANGL 

= 

AT AN (SLOPT(KK) ) 


ANGLl 

S 

ATANIYSU ,KK)/XS(1 ,KK) ) 


ANGL2 

= 

ATAN( YS(N,KK)/XS(N,KK) ) 


ANGLl 

= 

ANGL • , 5 * (ANGLl -TRAIL(KK)) 


ANGL2 


ANGL " . 5+ ( ANGL2 +TRAIL(KK)) 


T 1 

= 

TAN (ANGLl ) 


T2 


T A N ( A N G L 2 ) 


KK)), 

KK))) 
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CALL SPLIF C1,N,XP,YP,D1,D2,D3,1,T1,1,T2,0,0. 
CALL INTPL CI1,I2,X,Y,1,N,XP,YP,D1,D2,D3,0) 


-XI) 


XI 

» 

. 25*XS C 1 , KK ) 

A 

a 

SL0PTCKK)+CXSC1,KK) 

B 

= 

1. 

/CXSC1,KK) -XI) 

ANGL 

s 

PI 

+ PI 

U 

= 

1. 


V 

a 

0. 


M 

s 

11 

-1 

DO 52 

1 = 2 , M 




•5*SCAL*XCI)**2 
B* (XX -XI) 

YS(1,KK) +A*ALOG(I))/D 
SORT ( XX**2 +YY**2) 

ANGL + ATAN2 C CU^YY -V*XX),CU*XX 

XX 

YY 

SORT ( ( R +R)/SCAL) 

R*SIN C • 5*ANGL) 

SLOPT CKK)*CXSCN,KK) -XI) 
l./(XS(N,KK) -XI) 

0. 

1 . 

0. 

12 


+ 1 


XX 
D 

YY 
R 

ANGL 
U 

V 
R 

52 Y ( X ) 

A 
B 

ANGL 
D 

V 

M = 

DO 54 I=M,NX 
XX = 

D = 

YY = 

R = 

ANGL = 

U = 

V = 

R = 

54 Y ( I ) = 

Q = 

DO 62 1=2, NX 
62 SO ( I , K ) = SO ( 1 , K ) 

IF CKK.E0.K2) GO TO 71 
KK = K2 

P = R2 

GO TO 41 

71 DO 72 1=11,12 

72 I V ( I , K ) = 2 

M =11-1 

DO 74 1=2, M 
ZZ = COCK) 

IF (ZZ.GE.C0CKTE1) ) IV(I,K) = IVO 
74 CONTINUE 

M =12+1 

DO 76 I=M , NX 
ZZ = COCK) 

IF CZZ.GE.C0CKTE1)) IVCI,K) = IVO 
76 CONTINUE 

K2 = K2 -1 


,5*SCAL*XCI)**2 
B* C XX -XI) 

YSCN,KK) +A»ALOGCD)/D 
SORT C XX**2 + Y Y * * 2 ) 

ANGL +ATAN2C CU*YY -V*XX),CU*XX 

XX 

YY 

SURTCCR +R)/SCAL) 

R*SINC.5*ANGL) 

P/CG*CC) 


+Q*Y Cl) 


, IND) 


+ V*YY ) ) 


+ V*YY) ) 
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K = K +1 

IF (K.DE.KTE2) GO TO 21 
K l =2 

K2 = NZ 

81 DO 82 1=2, NX 

ZZ = COCK) 

IF (ZZ.bE.ZS(NC) . AND . ZZ . GE . CO C KTE1 ) ) IV(I,K) = IVO 

82 CONTINUE 

K = K +1 

IF CK.bE.K2) GO TO 81 
DO 102 K=K1,K2 
DO 104 1=2, NX 

IF CIVCI,K) .GT.O) GO TO 104 

IF (IV(I+1,K+1). GT.O. OR. IV(I-1,K+1). GT.O) IV(I,K) = IV1 
IF CIVCI+lrK-1). GT.O. OR. IVCI-1, K-l). GT.O) IV(I,K) = IV1 
104 CONTINUE 

102 IF CS0(LX,K) .LT.l.E-05) IVCLX,K) r 0 
DO 13 K=2,NZ 
DO 13 1=2, NX 

SI = SO (1+1 , K ) -50(1-1, K) 

SK = 50(1, K + l) -SO C I , K- 1 ) 

sx c i , k ) = Aicn* si 

SZ C I , K ) = C 1 C K ) * SK 

SXXCI,K) s (SO(I+1,K)-2.*SO(I,K)+SO(I-1,K)+A3(I)*SI)*A2(I) 
SZZ C I , K ) = (S0(1,K+1)-2.*S0(I,K)+S0(I,K-1)+C3(K)*SK)*C2(K) 

13 SXZ ( I , K ) = (S0(I+1, K+l) -S 0(1-1, K+l) -S 0(1+1, K-l)+S 0(1-1, K-l)) 

, ♦A1(1)*C1(K)40X*DX 

RETURN 
END 


♦DECK SINGPT 

SUBROUTINE SI NGPT ( X , Y , Nb , N , XSING , YSING , TRb , SLT) 

DIMENSION X(1),Y(1) 

RAD = 57.29577951308232 

NP x NL+1 

NM = NL-1 

CAbb XYSING(XCNb) , Y(Nb) ,X(NP) ,Y(NP) ,X(NM) ,Y(NM) , XSING, YSING) 
SbOPU = (Y(N)-Y(N-1))/(X(N) -X (N-l ) ) 

SbOPb = CY(1)-YC2))/(X(1)-X(2)) 

Sb r = .5+ csbopu +sbopb) 

THETAU a ATAN 2 C SLOPU , 1 . ) ♦RAD 

THETAb a AT AN 2 (SbOPb, 1 .) ♦RAD 

TRb = THETAb -THETAU 

RETURN 
END 
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♦DECK XYSING 

SUBROUTINE XYSING ( X 1 , Y 1 , X2 , Y 2 , X3 , Y 3 , XSI nG , YSING ) 

C FITS CIRCLE TO 3 POINTS NEAR LEADING EDGE AND FIND THE CENTER 

YA = ( Y2 + Y1)*.5E0 
XA = ( X2 + XI ) ♦ . 5EO 
YB = ( Y3 + Y 1 ) ♦ . 5E0 
XB a ( X3 + X 1 ) ♦ . 5£Q 
SL1 = -CX2 - XI) / CY2 - Yl) 

SL2 = -U3 - XI) / CY3 - Yl) 

XSING2 = (SL1 ♦ XA - SL2 ♦ XB + YB - YA) / (SL1 - SL2 ) 

XSING a (XSING2 + X1)*.5E0 
YSING2 a SL1 ♦ (XSING2 -XA) + YA 
YSING = (YSING2 ♦ Y1)4.5E0 
RETURN 
END 


♦DECK SPLIF 

SUBROUTINE SPLIF ( M , N , S , F , FP , FPP , FPPP , KM , VM , KN , VN , MODE , FQM , IND) 
C SPLINE FIT - JAMESON 

C INTEGRAL PLACED IN FPPP IF MODE GREATER THAN 0 

C IND SET TO ZERO IF DATA ILLEGAL 

DIMENSION SCI) ,F(1),FP(1),FPPC1),FPPP(1) 

IND = 0 

K = IABSCN -M) 

IF (K -1) 81,81,1 
IK s (N -M)/K 

I = M 

J = M +K 

DS = SCJ) -SCI) 

D = DS 

IF CDS) 11,81,11 

11 DF a CFCJ) -FCI))/DS 

IF CKM -2) 12,13,14 

12 U s .5 

V a 3 • ♦ C DF -VM ) /DS 

GO TO 25 

13 U =0. 

V s VM 
GO TO 25 

14 U = -1. 

V = -DS4VM 

GO TO 25 

21 I = J 

J = J + K 

DS = SCJ) -SCI) 

IF CD*DS) 81,81,23 
23 DF = CFCJ) -F Cl)) /DS 

B = l./CDS +DS +U ) 

U a B4DS 

V = &♦ C 6 . *DF -V) 
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non 


25 

31 

32 

33 

34 

35 

41 


51 

61 

71 


81 


♦ DECK 


FP ( I ) 

= U 

FPP(I) 

= V 

U 

= ( 2 . -U ) *OS 

V 

= 6 . ♦DF + DS*V 

IF (J - 

N ) 21,31,21 

IF (KN 

-2) 32,33,34 

V 

GO TO 35 

= ( 6 . ♦ V N -V)/U 

V 

GO TO 35 

= VN 

V 

= ( DS* VN + FPP ( I ) ) / ( 1 

B 

= V 

D 

= DS 

DS 

= SCJ) -SCI) 

U 

= FPPCI) -FP ( I ) * V 

FPPPCI) 

= (V -U ) /DS 

FPPCI) 

= U 

FP ( I ) 

= CF(J) -F Cl)) /DS - 

V 

= U 

J 

= I 

I 

= I -K 

IF CJ - 

M ) 41,51,41 

I 

= N -K 

FPPP(N) 

= FPPPCI) 

FPP ( N ) 

= B 

FP(N) 

= DF +D* ( FPP ( I ) +B 

IND 

= 1 

IF (MODE) 81,81,61 

FPPPCJ) 

= FQM 


. + FPC I ) ) 


DS*(V +U 


+B)/6. 


V 
1 
J 

DS 

U 

FPPPCJ) 

V 

IF (J 

RETURN 

END 



FPP C J ) 

J 

J +K 

SCJ) -S(I) 

FPP(J) 

FPPPCI) + . 5*DS* ( F ( 1 ) 
U 

71,81,71 


+ F ( J ) 


+U)/6. 


-DS*DS*(U 


♦ V)/12. ) 


INTPL 

SUBROUTINE INTPL (MI ,NI , S I , F 1 , M , N , 5 , F , FP , FPP , FPPP , MODE ) 
INTERPOLATION USING TAYLOR SERIES - JAMESON 
ADDS CORRECTION FOR PIECEWISE CONSTANT FOURTH DERI V I ATIVE 
IF MODE GREATER THAN 0 

DIMENSION SI(1) ,FIC1) ,S(1) ,F(1) ,FP(1) ,FPP(1) , FPPPCI) 

K n I ABS ( N -M) 

K = (N -M)/K 



NIN 

= N I 


D 

= S ( N ) -S(M) 


IF (D*CSI(NI) -SI(Ml))) 11,13,13 

11 

MIN 

= N I 


NIN 

= MI 

13 

KI 

= I ABS (NIN -MIN) 


IF (KI) 

21,21,15 

15 

KI 

= (NIN -MIN) /KI 

21 

II 

= MIN -KI 


C 

= 0. 


IF (MODE) 31,31,23 

23 

C 

= 1. 

31 

II 

= II +KI 


SS 

= SI(I1) 

33 

I 

= I +K 


IF (I 

-N) 35,37,35 

35 

IF CD* CSC I ) -SS)) 33,33,37 

37 

J 

= I 


I 

= 1 -K 


SS 

B SS -S(I) 


FPPPP 

= CnFPPP(U) -FPPP(I) )/(S(J) 


FF 

= FPPP(I) +. 25*SS*FPPPP 


FF 

= FPP(I) +SS*FF/3 , 


FF 

= FP ( I ) +.5^ssn-F 


FI ( I I ) 

= F ( I ) +SS^FF 


IF (II 

-NIN) 31,41,31 

41 

RETURN 



END 


♦DECK CPLOT 

SUBROUTINE CPLOT (I1,12,FMACH,X,Y,SU,SV,SW,SM,CP) 

C PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 

DIMENSION KODEC3),LINE(75),X(l),Y(l),SU(l),SV(l),Sw(l), 

. SM ( 1 ) , CPC 1 ) 

DATA KODE/1H ,1H+,1H^/ 

IWRIT = 6 

WRITE (IWRIT, 2) 

2 FORMAT (50HUPLOT OF CP AT EQUAL INTERVALS IN THE MAPPED PLANE/ 

1 8H0 X , 8H Y ,8H SU , 

2 8H SV ,8H SW , 8H SM ,8H CP ) 

FMACH2 = FM ACH*FMACH 

AAO = ( 1 . + . 2 *FM ACH2 ) 

CPO = (AA0*+3.5 -l.)/(.7*FMACH2) 

AAC = (1,+.2*FMACH2)/1«2 

CPC = ( AAC*^3 , 5 -l.)/(.7*FMACH2) 

DO 12 1=1,75 
12 LINE(I) = KODE(l) 

DO 22 1=11,12 

KC = 20. ♦(CPO -CPC) +20. 

KC = MAX 0(1 , KC) 
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KC = MIN0(75,KC) 

KK = 20 . * (CPO -CPU)) +20. 

KK = MAXO(t,KK) 

KK = M1N0 ( 7 5 , KK ) 

LINE(KC) = KODE ( 3 ) 

LINE(KK) = KODEC2) 

WRITECIWRIT , 6 1 0 ) X ( I) , Y (I) , SU ( I ) , SV ( I ) , SW C X ) , 
. SM(I) , CPU ), LINE 

LINE(KC) = KODE(l) 

22 LINE(KK) = KUOE(l) 

RETURN 

610 FORMAT ( 1 H ,7F8.3,75A1) 

END 


♦DECK FORCF 

SUBROUTINE FORCF ( I 1 , 12 , X , Y , CP , AL , CHORD , XM , CL , CD , CM ) 
C CALCULATES SECTION FORCE COEFFICIENTS 


DIMENSION 


RAD 

X 

ALPHA 

x 

CL 

= 

CD 

X 

CM 

X 

N 

X 

DO 12 1=11 

# N 

DX 

X 

DY 

r 

X A 

X 

YA 

X 

CPA 

X 

DCL 

X 

DCD 

X 

CL 

X 

CD 

X 

12 CM 

X 

DCL 

s 


X(1),Y(1),CP(1) 

57.2967795130823 

AL/RAD 

0. 

0. 

0. 

12 -1 

(XCI+1) -X ( I ) ) /CHORD 

(Yu+i) •¥(!)) /chord 

(,5»(X(I + 1) +X C 1 ) ) -XM ) /CHORD 
• b*U(I + l) + Y (I ) ) /CHORD 
. 5* CCP U + l ) +CPID) 

-CPA ♦DX 
CPA^DY 
CL +OCL 
CD +DCD 

CM +DCD* Y A -DCL^XA 
CL^COS( ALPHA) -CD^SIN ( ALPHA ) 


CD = CL+SIN ( ALPHA) +CD^C05( ALPHA) 

CL = DCL 

RETURN 

END 


♦DECK TOTFOR 

SUBROUTINE TOTFOR l KTE 1 , KTE2 , CHORD , SCL , SCI) , SCM , CO , XC , 
1 CL,CD,CMP,CMR,CMY) 

C CALCULATES TOTAL FORCE COEFFICIENTS 
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DIMENSION 


CHORD ( l ) , SCL Cl), SCD ( 1 ) , SCM ( 1),C0(1),XC(1) 

SPAN 

3 

C0CKTE2) -C0CKTE1) 

CL 

= 

0. 

CD 

s 

0. 

CMP 

3 

0. 

CMR 

3 

0. 

CMY 

*+ 

0. 

S 

3 

0. 

N 

3 

KTE2 -1 

DO 12 KsKTEl 

, N 

DZ 

3 

. 5* CCO C K + l ) -COCK)) 

AZ 

3 

.5MC0CK + 1) +C0CK)) 

CL 

3 

CL +DZ^CSCLCK+l)»CHORDCK+l) +SCL C K ) ♦CHORD C K ) ) 

CD 

3 

CD +DZtCSCDCK+lMCHURDCK+l) +SCD C K ) ♦CHORD C K ) ) 

CMP 


CMP +DZ^ C CHORD C K+l )^CSCMCK+1) ♦CHORD C K+l) 

1 


-SCLCK + imCCK + 1) ) 

2 


+CHORDCK)^CSCMCK)^CHORDCK) 

3 


-SCL C K ) *XC C K ) ) ) 

CMR 

3 

CMR +AZ*DZ*CSCLCK+l)^CHORDCK+l) +SCL C K ) ♦CHORD ( K ) ) 

CMY 

3 

CMY +AZ^DZMSCDCK + l)+CHORDCK+l) +SCD C K ) ♦CHORD C K ) ) 

12 S 

3 

S +DZ*CCRQRDCK+1) +CHORDCK) ) 

CL 

s 

CL/S 

CD 

3 

CD/S 

CMP 

3 

CMP^SPAN/S^2 

CMR 

3 

C CMR +CMR)/CS + SPAN) 

CMY 

RETURN 

END 


CCMY +CMY)/ CS^SPAN) 


♦ DECK 

C 

C 

♦ CALL 


1 
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PSURE 

SUBROUTINE PSURE ( I PLOT ,K,X,Y,CP,I1,I2,CL,CD,CM) 

GENERATES PLOT FOR PRESSURE DISTRIBUTION OVER BLADE SECTION 
AT EQUAL INTERVALS IN THE MAPPED PLANE 


BLANK 

DIMENSION R(IOO) ,D1 (ISO) ,D2(150) ,D3( 150) 
DIMENSION XC1),YC1),CPC1) 

IF (IPLOT) 1,11,101 
CALL VERSA ( 20 ) 

CALL BGNPL(-l) 

IPLOT = 0 

CALL PHYSOR C 0 . , 0 . ) 

CALL TITLE ( 1 H ,0,1H ,0,1H ,0,8. ,10.5) 


CALL GRAPH ( 0 
ZSO = 

ZS = 

VROT = 

VTAN = 

SMACH(K) = 

T 1 = 

PSID = 


«, 1 ., 0 ., 1 .) 

CCO(K) +HINGE)/(C0(KTE2)+HINGE) 

COCK) + HINGE 

OMEGA*ZS 

VROT +FMACH»SIN(PSI)+COSCALPHA) 
VTAN 

l./(.7+SMACH(K)*^2) 

PSI ♦R AD 
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ENCODE ( 45 , 4 , R ) PSID, FMACH,TMACH 
4 FORMAT ( bHPSI s , f 7 . 1 , 3X , bHFMACHs , F7 . 4 , 3 X , 6HTM ACH= , F7 . 4 ) 
CALL MESSAG(R, 45,1. 5,1.) 

ENCODE (45 , 1 5 , R ) ZSO , SMACH ( K ) , TILT 

15 FORMATCbHZS = , F7 . 4 , 3X , 6HSMACH= , F 7 . 4 , 3X , bHAL =,F7.4) 
CALL MESSAG(R,45,1.5,0.75) 

ENCODE (45,16,R) CL, CD, CM 

16 FORMAT ( 6HCL = , F 7 . 4 , 3X , 6HCD « , F7 . 4 , 3X , 6HCM =,F7.4) 

CALL MESSAG(R,45, 1 .5, .5) 

ENCODE(2, 17 r R) 

17 FORMATC2HCP) 

CALL MESSAGCR, 2, 1 .4,5.25) 

C DRAW AIRFOIL 

XMAX = X(I1) 

XMIN = X(ll) 

DO 22 I =11,12 

XMAX = AMAXl(XCI) ,XMAX) 

22 XMIN = AMINKX(I) ,XMIN) 

SCALE = 5. /(XMAX -XMIN) 

XOR = 2. 

YOR = 2. 

N s 12-11+1 

DO 24 J=1,N 

D 1 ( J ) = SCALE*(X(J+I1— 1)-XMIN) +XOR 

24 D2(J) = SCALE* Y (J+Il-1) +YOR 

CALL CUHVECDl ,D2,N,0) 

CPMAX = 0. 

IMAX = XI 

DO 25 1= 11,12 
ABSCP = CP(I) 

IF(ABSCP.LE.CPKAX) GO TO 25 
CPMAX = AbSCP 

IMAX = 1 

25 CONTINUE 

YOR = YOR + 3. 

C CPC IS CRITICAL PRESSURE COEFFICIENT 

AAC = (l.+.2*SMACH(K)**2)/1.2 

CPC = ( A AC** 3 . 5 - 1. )*T1 

IF(ABS(CPC) .GT.1.2) GO TO 50 
CPCM = Y0K-2.54CPC 

CALL STRTPT(2. ,CPCM) 

CALL C0NNPTC3. ,CPCM) 

50 N a IMAX -II +1 

DO 32 J=1,N 
D3(J) = D 1 ( J ) 

32 D2 ( J ) = YOR -2.5*CP(J+I1-1) 

CALL MARKER ( 4 ) 

CALL CIJRVE(D3,D2,N,0) 

N = 12 - IMAX + 1 

DO 34 J= 1 , N 

D3 ( J ) = DHJ + IMAX-I1) 

34 02 (J) = YOR-2.5*CP( J+1MAX-1 ) 

CALL MARKER ( 3 ) 
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CALL CURVE(D3,D2,N,0) 

CALL ENDGR(O) 

C DRAW CP AXIS 

CALL OREL ( 2 • , 2 . ) 

CALL TITLEC1H ,0,1H ,0,1H ,1,6. ,6.) 
CALL YAXANGCO.) 

CALL GRAPH(0.,1. ,1.2,-.4) 

CALL ENDPL(O) 

101 RETURN 
END 


♦DECK THREED 

SUBROUTINE 'THREED ( I PLOT , SU , SV , SW f SM , CP , X , Y , TYTLE , CHORD , 

1 CL, CD,CHORDO,SCL,SCD,SCM) 

C GENERATES PLOT FOR PRESURE DISTRIBUTIONS OVER BLADE 

♦CALL BLANK 
♦CALL A 

DIMENSION X(1),Y(1),SU(1),SV(1),SW(1),SM(1),CP(1), 

. SCL ( 1 ) ,SCD (1),SCM(1), CHORD ( 1 ), TYTLE C 1 ) , 

. XD(200) , YD (200) ,CPD(200) ,R(80) 

IF ( I PLOT )1,11,101 
1 CALL VERSA C 20 ) 

CALL BGNPL(-l) 

IPLOT = 0 

11 CALL PHYSOKCO. ,0. ) 

CALL TITLE(1H ,0,1H ,0,1H ,0,8. ,10.5) 

CALL GRAPH(0.,1.,0.,1.) 

SPAN s C0CKTE2) -CO(KTEl) 

AR = SPAN/CHORDO 

SCALXX * 2.5/CH0RD0 

SCALZZ = 5. /SPAN 

SCALPP = -1.25 

TX = 3.5 

XOR = 4.5 -SCALXX ♦XC(KTEl) 

YOR = 3.75 

DO 6 K= KTE1 , KTE2 

11 = ITEl(K) 

12 = ITE2CK) 

CALL VELO(K,SU,SV,SW,SM,CP,X,Y) 

CHORD ( K) b XC11) -X(LX) 

CALL FORCF ( 1 1 , I 2 , X , Y , CP , TILT, ChORD ( K ) , XC ( K ) , SCL (K ) , SCD ( K ) , SCM ( K ) ) 
SY = SCALZZ* ( CO ( K ) -CO(KTEl)) +YOR 

DO 7 1= II, LX 
JsI-Il+1 

XD(J) r SCALXX*X ( I ) +XOR 

7 CPD(J) s SCALPP*CP ( I ) +SY 

N = LX -II +1 

CALL CURVE(XD,CPD,N,0) 

DO 8 I=LX ,12 
0 = I-LX+1 
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XD(J) = SCALXX*X ( I ) + XOR -XX 

8 CPD(J) = SCALPP*CP ( I ) +SY 

N =12 -LX +1 

CALL CURVE(XD,CPD,N, 0) 

6 CONTINUE 

CALL MESSAG(49HUPPER SURFACE PRESSURE LOWER SURFACE PRESSURE, 

. 49,1.5,1.5) 

CALL MESSAGCTYILE, 100,1. ,1.) 

ENCODE l 45, 3, R) F MACH, T MACH, TILT 

3 FORMAT(6HFMACH=,F7.4,3X,bHTMACH=,F7.4,3X,6HALPHA=,F7.4) 

CALL M ES S AG (R, 45,1. ,0.75) 

CALL TOTFOKCKTEl ,KTE2, CHORD, SCL,SCD,SCM, CO, XC, 

1 CL,CD1 ,CMP,CMR,CMY) 

CD = CD 1 

PSIO = PSI ♦RAD 

ENCODE ( 45 , 4 , R ) PSID,CL,CD 

4 FORMAT ( 6HPSI =, F7 . 1 , 3X , 6HCL = , F7 . 4 , 3X , OHCD =,F7.4) 

CALL MESSAGCR,45,1. ,0.5) 

CALL ENDPL(O) 

101 RETURN 
END 


♦DECK ROTORtt 

SUBROUTINE ROTORB ( I PLOT , SU , SV , SW , SM , CP , X , Y ) 

C GENERATES PLOT FOR ROTOR BLADE GEOMETRY 

♦CALL BLANK 
♦CALL A 

DIMENSION X(l), Y(1),SU(1),SV(1),SW(1),SM(1),CP(1), 
. D 1 C 2 0 0 ) ,D2(200) ,D3(200) ,D4(200) ,D5(200) , 

. XSMAXC50) ,XSMIN(50) ,ZSTAT(50) ,R(80) 

IF ( IPLOT) 1,11,101 
1 CALL VERSA ( 20 ) 

CALL BGNPL(-l) 

IPLOT = 0 

11 CALL PHYSOR ( 0 . , 0 . ) 

CALL TITLEC1H ,0,1H ,0,1H ,0,8. ,10.5) 

CALL GRAPHC0.,1.,0.,1.) 

SPAN = C0CKTE2) -CO(KTEl) 

SCALZZ = 7. /SPAN 

CALL MESS AG ( 2 1 HONER A SWEPT TIP BLADE , 2 1 , 2 . , 1 . ) 

DO 6 K= KTE 1 , KTE2 

11 x lTEl(K) 

12 = ITE2CK) 

CALL VELO(K,SU,SV,SW,SM,CP,X,Y) 

XMAX = -10. 

XM IN = 10. 

DO 7 1=11,12 

XMAX = AMAX1 (XMAX,X(I) ) 

7 XMIN = AM1N1(XMIN,X(I)) 

IFCK.EQ.KTE1 ) XSTAT= XMIN 
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o n 


XSMAX(K) = SCALZZ*(XMAX -XSTAT) 

XSMIN(K) = SCALZZ* ( XM1N -XSTAT) 

ZS = CO ( K ) -CO(KTEl) 

ZSTAT(K) = SCALZZ*ZS 
6 CONTINUE 

DO 8 K=KTE1 , KTE2 
KK = K-KTE1+1 

Dl(KK) = 2.+XSMAXCK) 

D2 (KK) s 2.+XSMINCK) 

8 D3 ( KK ) = 2.+ZSTATCK) 

N = KTE2-KTE1+1 

CALL CURVECDl ,D3*N,0) 

CALL CURVE(D2,D3,N,0) 

C DRAW AIRFOIL 

DO 21 KK=KTE1,KT£2 
K = KK-KIEl +1 

11 = 1TE1(KK) 

12 = ITE2 ( KK ) 

CALL VELO(KK,SU,SV,SW,SM,CP,X,Y) 

N = I2-I1+1 

DO 24 J=1 ,N 

D4( J) = SCALZZMX(J + 11-1 )-XSTAT) +2. 

24 D5( J) = SCALZZ*Y(J+I1-1 ) +D3CK) 

CALL CURVE(D4,D5,N,0) 

21 CONTINUE 

CALL ENDPL(O) 

101 RETURN 
END 


♦DECK VELO 

SUBROUTINE VELG(K,SU r SV,SW,SM,CP,X,Y) 
CALCULATES SURFACE VELOCITY 
CP SCALED BY FAR FIELD SOUND SPEED 
♦CALL BLANK 
♦CALL A 


DIMENSION 

AAO 

= 

SU(1) ,SV(1) ,SW(1),SM(1) ,CP(1),X(1) ,Y(1) 

i. 

11 

— 

ITEl(K) 

12 

x 

ITE2CK) 

ZS 

X 

COCK) + HINGE 

VROT 

X 

OMEGA*ZS 

VTAN 

X 

VROT +FMACH*SIN(PSI)*COS( ALPHA) 

SMACH(K) 

X 

VTAN 

T 1 

X 

1»/(*7*SMACH(K)**2) 

DO 12 1=11 

f 

12 

XI 

X 

AO ( I ) 

Y 1 

X 

SO ( I , K ) 

X1X1 

X 

X1*X1 

Y1Y1 

X 

Y 1 ♦ Y 1 

HH 

X 

X1X1 + Y1Y1 
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DHH 

= 

1./ HH 

XB 


. 5* ( X 1 X 1 -Y1Y1) 

YB 

= 

X 1 * Y 1 

XS 

= 

XC(K) +XB*SCAL 

X 1 XB 

s 

XI *DHH 

XI YB 

= 

Yl» DHH 

X 1 ZB 


-XZ(K) * X 1 X B -YZ(K) *X1YB 

Y 1 ZB 

= 

XZ ( K ) ♦ X 1 YB - YZCK) *X1XB 

YXB 

— 

“(XI YB +XlXB*SX(I f K)) 

YYB 

— 

XlXb -X1YB*SX(I,K) 

YZ8 

s 

Y 1 ZB -X 1 ZB*SX ( I , K ) -SZ(I,K) 

GI 

= 

G(I+1,2,K) -G(I-1,2,K) 

GJ 

s 

2 . * ( G ( I , i , K ) -G ( I , 2 , K ) ) 

GK 

= 

G(I,2,K+1) -G(I,2,K-1) 

GX 

s 

A1 ( 1 ) *GI 

GY 

= 

B1(2)»GJ 

GZ 

= 

Cl (K) *GK 

U 

= 

(GX+X1XB+GY*YXB)/SCAL 

V 

= 

(GX*X1 YB+GY*YYB)/SCAL 

W 

r 

(GX*XlZB+GY»YZBtGZ)/SCAL 

Q Q 

s 

U*U+ V*V + W*W 

UF 

= 

DMEGA*ZS +CAC 

VF 

s 

SA 

WF 

= 

“ (UMEGA*XS +CAS ) 

TERMS 

ss 

U*UF + V*VF + W*WF 

FIT 

= 

TERMS 

AA 


DIM(AAO,.2»QO+.4*FIT) 

UB 

= 

U +UF 

VB 

= 

V + VF 

MB 


W +WF 

UUB 

= 

UB^UB 

V V B 

= 

VB* VB 

W MB 

s 

MB*wB 

QQK 

s 

UUB +VVB +WMB 

SU(1) 

= 

UB 

SVCl) 

= 

VB 

5W ( I ) 

= 

MB 

SM ( I ) 

=: 

SGRTCOUR/AA) 

CP(I) 

s 

(AA**3.5-1.)*X1 

X(I) 

= 

XS 

12 YCI) 

r 

YC ( K ) +SCAL* YB 

RETURN 



END 




♦ DECK 

ESTIM 


SUBROUTINE ESTIM 

C 

INITIALIZATION FOR STEADY CALCULATION 

♦ CALL 

BLANK 

♦ CALL 

A 


LX a NX/2 +1 
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c 


DX 


= 2. 

/FLOAT ( NX) 

DSUM 


= 1. 

/ ( A2 ( LX ) +B2 ( 2 ) ) 

WATY 


a B2(2)*DSUM 

ini ATX 


a A2 (LX ) ♦DSUM 

AAO 


a 1. 


MX 


a NX 

+ 1 

MY 


a NY 

+ 2 

MZ 


a NZ 

+ 1 

DO 17 

Ja 

1 ,MY 


DO 17 

Ka 

1,MZ 


DO 17 

la 

1,MX 


GCI, J 

,K) 

a 0. 


17 CONTINUE 




SURFACE CONDITION 


DO 23 K=2,NZ 

IF(ITE2(K) .EQ.MX) GO TO 23 
ZS = CO ( K ) +HINGE 

1X1 = ITE1(K) 

1X2 s ITE2 (ft ) 

DO 22 1=1X1, 1X2 
XI = AO ( I ) 

Y 1 = S0(I,K) 

X1X1 = XI ♦ Xl 

Y1Y1 = Y1 ♦*! 


HH 

DHH 

Xb 

XS 

X 1 XB 
XI YB 
X 1 ZB 

Y 1 ZB 
GI 
GK 
GX 
GZ 
UF 
VF 
WF 
YXB 

Y YB 
YZB 
X1YS 
YYS 
RHS 

G(I,1,K) 

22 CONTINUE 

23 CONTINUE 
RETURN 
END 


= X1X1 +Y1Y1 
= l./HH 

a .5*(X1X1 -Y1Y1) 
a XC ( K ) +XB*SC AL 
= XI ♦ DHH 
= Y 1 ♦DHH 

a -XZCK) ♦ X 1 X B -YZ(K) ♦ XIYB 
= XZ(K) ♦ XIYB -YZCK) + X1XB 
= G(I+1,2,K) -G(I-1,2,K) 

= G(I,2,K+1) -G(I,2,K-1) 

= A1 C I ) *GI 
= C 1 ( K ) ♦GK 
a OMEGA^ZS +CAC 
= SA 

= -COMEGA^XS +CAS ) 

= -(XIYB +X 1 XB^SX ( I , K ) ) 

= X1XB -X1YB*SX(1,K) 

a Y 1 ZB -X1ZB*SX(I,K) -SZ(I,K) 

a XIXB^YXB +X1YB4YYB +X1ZB*YZB 

= YXB*YXB +YYB*YYB +YZB^YZB 

B ( U F ♦ Y X B +VF*Y YB + WF»Y ZB) ♦SCAL 

= G(I,3,K) + (RHS +X1YS+GX + Y2B«GZ) / ( YYS^Bl ( 2 ) ) 
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♦DECK REFIN 

SUBROUTINE REFIN 
♦CALL BLANK 
♦CALL A 

LX = NX/2 +1 

A AO = 1. 

DSUM = 1./CA2CLX) +b2 C 2 ) ) 

WATY = B2 ( 2 ) ♦ DSUM 

WATX a A2(LX)*OSUM 

DX = 2. /NX 

MX = NX +1 

MY = NY +2 

MZ = NZ +1 

M X 0 = NX/2 +1 

MYO = NY/2 +2 

MZO = NZ/2 +1 

DU 1 MK=l,MZO 
K =MZO +1 -MK 

KK = CK-1 )^2 +1 

DO 1 M,J = 2,KYO 
J =MYO +2 -MJ 

JJ = (J-2)*2 +2 

DO 1 MI=i,MX0 

I a MXO +1 “MX 

II = C I “ 1 ) ♦ 2 +1 

1 G(U,JJ,KK)=G(I,J,K) 

DC) 2 K=1,MZ,2 

DO 3 J=2,MY,2 
DO 3 1=2, MX, 2 

3 G(i,J,K) = , 5*(G(I+1,J,K) +G ( 1-1, J,K)) 
DO 4 1=1 , MX 

DO 4 J=3,MY,2 

4 G(I,J,K) = .5^(G(I, J.+ 1 , K ) +G( I , J-l , K) ) 

2 CONTINUE 

DO 5 K=2 , MZ , 2 
DO 5 J=2 , M Y 
DO 5 1=1, MX 

5 G(I,J,K) = . 5*(G(I,J,K+1) +G(I,J,K-D) 

DO b K=2 , NZ 

1X1 = ITEl(K) 

1X2 = ITE2CK) 

IF C 1X2 . EQ . MX ) GO TO 7 
ZS = CO(K) +H I NGE 

C WING CONDITION 

DO 10 1=1X1, 1X2 
XI = AO(I) 

Yt = SOCI ,K) 

X1X1 = XI ♦ X 1 

Y1Y1 = Y 1 ♦ Y 1 

HH = X1X1 +Y1Y1 

DHH = 1. /HH 

XB = , &♦ ( X 1 XI -Y1Y1) 

XS = XCCK) +XB+SCAL 

X 1 XB = XI ♦DHH 

X 1 Yb = Y 1 ♦DHH 
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X 1 ZB 

= -XZ(K) ♦XIXB -YZCK) ♦XIYB 


Y 1 ZB 

= XZ(K) ♦XIYB -YZ(K) ♦XIXB 


GI 

= G(I+1,2,K) -G(I-1,2,K) 


GK 

= G ( 1 , 2 , K+ 1 ) -G(I,2,K-1) 


GX 

= A1 ( I ) ♦GI 


GZ 

= Cl (K ) ♦GK 


UF 

= OMEGAtZS +CAC 


VF 

= SA 


WF 

= -(OMEGA+XS +CAS) 


YXB 

= “(XIYB +X 1 XB^SX ( 1 , K ) ) 


YYB 

= XIXB -X 1 YB*SX ( I , K ) 


YZ8 

= Y 1 ZB -X1ZB*SX ( 1 , K ) -SZ(I,K) 


X1YS 

= X1XB+YXB +X1 YB^YYB +X1ZB^YZB 


YYS 

= YXB+YXB +YYB^YYB +YZB^YZB 


RHS 

= (UF^YXB +VF+YYB +*F^YZB) ♦SCAL 

10 

G(I,1,K) 

= G(1,3,K) +CRH5 +X1 YS^GX + Y7.B+GZ) / ( YYS^Bl ( 2 ) ) 


E 

= G(IX2,2,K)-G(IX1 ,2,K) 


IX 

= 1X2 +1 


DO 8 I=IX 

,MX 


M 

= NX +2 -I 


G ( 1 , 1 , K ) 

= G(M, 3,K) +E 

8 

G(M,1,K) 
GO TO 6 

= G(1,3,K) -E 

7 

G ( LX , 2 , K ) 

= G(LX,3,K) + WATY+G(LX-1,2,KMWATX 


DO 9 I =LX 

f MX 


M 

= NX +2 -I 


G(I,2,K) 

= G(M,2,K) 


G(I,1,K) 

= G(M,3,K) 

9 

G(M, 1 r K) 

= G(1,3,K) 

6 

CONTINUE 

RETURN 

END 


♦ DECK 

RELAX 



SUBROUTINE RELAX 

♦ CALL 

BLANK 


♦ CALL 

A 


♦ CALL 

FLO 



DIMENSION 

C(131),D(131),GM(129,18,33), 
AB(129),AC(129),AA(129) ,GGR(129) ,R(129) , 
HH(129) , XXISQ29) , YYS (129) ,X1YS(129) , 

XI ZB (129) , YZBC129) ,GI(129) ,GJ(129) ,GK(129), 
G 11(129) , GJJ(129) ,GKK(129) , 
GIJ(129),G1K(129) ,GJK(129) , 
UUR(129),VVR(129),KhR(129), 
UVR(129),UWR(129),VKR(129), 

UK ( 1 29 ) ,VR(129),wR(129) 


* T 1 

= DX^DX 


01 

= 2. /PI 


Q2 

= 1./R2 


FR 

= 0. 


IR 

= o 


JR 

= 0 


KR 

= 0 


GD 

= o. 


IG 

= 0 
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70 

303 


C 


31 


C 

403 


JG =0 

KG =0 

NS = 0 

C(l) = 0. 

D ( 1 ) =0. 

DO 70 K=1,MZ 
DO 70 J=1,MY 
DO 70 1=1, MX 
GM ( I , J , K ) = G(I,J,K) 

CONTINUE 

DO 103 K = 2 , NZ 

ZS = CO ( K ) +HINGE 

FOR FIXED WING FLOW J=NY+1 

J = NY 

13 = NX 

BC = Tl*81(J)*Cl(K) 

INTERIOR 

DO 400 1= 2,13 

AB(I) = T1*A1 (I) *B1(J) 

AC C 1 ) = T1 * A 1 ( I ) * C 1 C K ) 

XI = AO ( I ) 

Y 1 = BO ( J ) +S0 ( I , K ) 

X1X1 = XI *X1 

Y1Y1 B Y 1 * Y 1 

HH(I) = X1X1 +Y1Y1 

DHH = 1. /HH ( I ) 

XB = .5*(X1X1 -Y1Y1) 

XS = XCCK) +XB*SCAL 

X 1 XB = XI ♦DHH 

X1YB = Y 1 *DHH 

X1ZB(I) = -XZ(K) *X1XB -YZIK) *X1YB 

Y 1 ZB = XZCK) ♦XIYB -YZ(K) *X1XB 

YXB = -CX1YB +X 1 XB*SX ( 1 , K ) ) 

YYB = X 1 XB -X1YB*SX(I,K) 

YZB(I) = Y 1 ZB -X1ZB(I)*SX(I,K) -SZ(I,K) 


GICI) 

GJCI) 

GK ( I ) 

GX 

GY 

GZ 

U 

V 

W 

00 

UF 

VF 

WF 

TERMS 

FIT 

AA(I) 

UB 

VB 

WB 

U1JB 


G(I+1,J,K) -G(I-1,J,K) 
GC1,J41,K) -G(I,J-1,K) 
GCI,J,K+1) -GCI,J,K-1) 

A 1 t I ) *G I ( I ) 

B1 ( J ) ♦G J ( I ) 

C1(K)*GK(I) 
(GX*X1XB+GY*YX6)/SCAL 
(GX*X1YB+GY*Y YB)/SCAL 
(GX*X1ZB(I)+GY*YZB(I)+GZ)/SCAL 
u*u + v*v +w*w 

OMEGA*ZS +CAC 
SA 

-(UMEGA*XS +CAS ) 

U*UF + V*VF + W*WF 
TERMS 

D I M ( A A 0 , . 2 * QQ + . 4 *F I T ) 

U + UF 
V + V F 
w + W f 
UB*UB 
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U\ & W NJ 


VVB 

WWB 

UVB 

UWB 

VWB 

UR ( I ) 

VRCI) 

WR ( I ) 

UURCI) 

VVR(I) 

WWRCI J 

UVR(I) 

UKR(I) 

VWRCIJ 

GGk(I) 

XXIS(I) 

XIYS(I) 

YYSCI) 

X1XBXB 

X1XBYB 

BCHI 

BPSI 

BLAMD A 

BSIGMA 

FA 

FB 

FC 

FD 

FE 

FF 

FA A 

FBB 

FCC 

FDD 

FEE 

RL 

1 

2 

RM 

1 

2 


RN 

400 R ( I ) 

DO 401 1= 

GIIU) 

GJJ(I) 

GKK(I) 

GIJ(I) 

GIKU) 

GJKCI ) 

401 CONTINUE 


= VB*VB 
= WB*WB 
= U B * V b 
= U B * tf B 
= VB*fcb 

= X1XB*U8 +X1YB*VB +X1ZB(I)*W6 
= YX6+U0 + Y Y B ♦ V B +YZB(I)*WB 
= WB 

= UR(I)*UR(I) 

= VR(I)*VR(I) 

= WWB 

= UR(I)*VR(I) 

= URCI)*WR(I) 

* VR(I)*WR(1) 

= UUB +VVB + WWB 

= X 1 XB*X 1 XB + X 1 YB*X 1 YB + X 1 ZB ( I ) *X 1 ZB ( I ) 

= X1XB*YXR +X1 YB*YYB +X 1 ZB ( I ) * YZB ( I ) 

= Y X b ♦ Y X 8 +YYB* Y YB +YZB ( 1 ) * YZB C I ) 

= -X1*(HH(I) -4.*Y1Y1)*DHH**3 
= Y1*(HH(1) -4.»X1X1)*DHH**3 

= XZ C K ) * X 1 XBXB +YZ(K)*X1XBYB 
= XZ ( K ) *X 1 XBYB -YZ(K)*X1X8X6 
= SX(I ,K)*X1XBXB +X1XBYB 
= SX(I,K)*X1XBYB -X1XBXB 
= XZ ( K ) * BCHI +YZ(K)*BPS1 
= YZIK)*BCHI -XZ ( K ) ♦BPS 1 
= XZ(K)*BLAKDA +YZCK) ♦BSIGBA 
= XZ ( K ) ♦BSIGM A “YZ( K ) *BLABDA 
= XZZ ( K ) *X 1 XB +YZZ(K ) *X1 YB 
= XZZ(K)*YXb +YZZ(K)*YYB 
= FA-F£*SCAL 
= OMEGA*(CAS“WFO*SCAL 
= OMEGA*(CAC+UF)*SCAL 
= FB-FA*SX(I ,K) 

= FDD-FF*SC AL 

= T1*(-FAA*(WWB-AA(I ) ) -X1XBXB*(UUB-VVB) 

+2 • * ( BCHI*UWB “X1XBYB*UVB +BPSI*VWB) 

+FB0*X 1 XB +FCC*X1ZB ( I ) ) 

= T1*(-FEE*(WW8-AA(1)) 

+BLARDA*(UUB-VVB)+2.*(UVB*BSIGMA-UWB*FC-VWB*FD) 

+ FBB* YXB +FCC* YZB ( 1 ) ) 

“AA(IJ*(XX1S(I)*SXX(1,K) + SZZ ( I , K ) 
+2.»X1ZB(I)*SXZ(1,K) ) 

♦ UUR(I)*SXX(I , K ) +WWB*SZZ(I,K) 

+2.*UR(I)*WB*SXZ(I,K) 

= Tl+FCC 

s RL*GX +RM*GY +RN+GZ 
2,13 

= G(I+l,J,K)-2.*GCI,J,K)+G(I-l,d,K) +A3 C I ) *GI ( I ) 

= G(I,J+1,K)-2.*G(I,J,K)+G(I,J-1,K) +B3 ( J ) *GJ ( I ) 

= G(J,J,K+l)-2.*G(I,d,K)+G(I,J,K-l) +C3(K)*GK(I) 

= G(I+1,J+1,K)-G(I+1,J-1,K)-G(I-1,J+1,K)+G(I-1,J-1,K) 
= G(I+1,J,K+1)“G(I+1,J,K“1) “G (I“1,J,K+1)+G(I“1,J,K“1) 
= G(1,J+1 ,K+1 )-G(I , J+l ,K-1)-G(1,J-1 ,K+1)+G(I,J-1,K-1) 
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DO 8 1=2 ,13 

SIGNX = SIGN ( 1 . ,Uk(I) ) 

SIGN* = SIG N(1.,VR(I)) 

SIGNZ = SIGN ( 1 . , WR( I ) ) 

AXT = BETA*UR(I)*A1(I) 

AYT = 6ETA*VR(1)*B1(J) 

AZT = BETA*WR(I)*C1(K) 

LL = IFIX(SIGNX) 

IM = I-LL 

IMM = IM-LL 

LL s IFIX(SIGNY) 

J M = J — LL 

JMM = JM-LL 

LL = lFIX(SIGNZ) 

KM = K-LL 

K MM = KM-LL 

IFCOQR(I) .GE.AA(I) ) GO TO 9 

AXX = A2(I)*(UUR(I)-AA(I)*XX1S(I)) 

AY Y = B2(J)*(VVR(I)-AA(I)*YYS(I)) 

AZZ = C2(K)*(WWR(I)-AACI)) 

AX Y = 2.*ABU)*(UVR(I)-AA(I)*X1YS(I)) 

AXZ = 2.*AC(I)*(UWR(I)-AA(I)*X1ZB(I)) 

AYZ = 2.*RC*(VwR(l)-AA(l)*YZB(I)) 

Y I = -(AXX*GII(1)+AYY*GJJ(I)+AZZ*GKK(I) 

. +AXY*GIJ(I)+AXZ»GIKCI)+AYZ^GJK(I))+R(I) 

Cl = -AXX 

B I = -AXX 

A I = AXX + AXX +U1*(AYY +AZZ) 

GO TO 10 

C TYPE DEPENDENT DIFFERENCING 

9 NS = NS +1 

BXX = A2C1) * (QOH(I) ♦XXIS(I) -UUR(D) 

BY Y = B2(J) * COORCI) ♦YYS(l) -VVR(I)) 

BZZ = C2CK) * (OQR(l) -WRRCI) ) 

BXY = 2.*AB(I) * (GQR(l) ♦XIYS(I) -UVR(I)) 

BXZ = 2.»AC(I) * (GOR(I) ♦XIZB(I) -UWR(I)) 

BYZ = 2.*t>C * (OQR(I) ♦ YZB ( I ) -VWR(I)) 

DELTA = BXX*GII(I)+BYY*GJJ(I)+BZZ*GKK(I) 

. +BXY*GIJ(I)+BXZ*GIK(I)+BYZ*GJKCI) 

IFCIMM.LT.l.OR.IMM.GT.MX) GO TO 11 
GII(I) = G(I,J,K)-2.*G(1M,J,K)+G(IMM,J,K) 

. + 2.»A3(I)+SiGNX*lG(I,J,K)-G(IM, J,K) ) 

11 IF(JMM.LT.l.OR.JMM.GT.MY) GO TO 12 

GJJ(I) = G(I, J,K)-2.*G(I,JM,K)+G(I,JMM,K) 

. + 2.*83(J)*SIGNY*(G(I,J,K)-G(I, JM,K)) 

12 IF(KMM.LT.l.OR.KMM.GT.MZ) GO TO 13 

GKK(I) = G(I, J,K)-2.^G(I,J,KM)+G(I,J,KMM) 

. + 2.*C3(K)*SIGNZ*(G(I, J , K ) -G ( I , J , KM ) ) 

13 GIJ(I) = G(I,J,K)- G(IM,J,K)- G(1,JM,K)+ G(IM,JM,K> 

GIKCIJ = G(I,J,K)- G(IM,J,K)- G(I,J,KM)+ G(IM,J,KM) 

GJK(IJ = G(1,J,K)- GCI,JM / K)- G(I,J,KH)+ G(I,OM,KM) 

AXX = (JUK ( I ) * A2 ( I ) 

A Y Y = VVRtl)*B2(J) 

AZZ = WWk(I)*C2CK) 

AX Y = 8. ♦SIGNX ♦SIGNY ♦AB(I) ♦UVR(I) 
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AXZ = 8. +SIGNX ♦SIGNZ ♦AC(I) *UWR(I) 

AYZ = 8. ♦SIGN Y ♦SIGNZ ♦ BC ♦VWR(I) 

GSS = AXX»GII(I)+AYY»GJJCI)+AZZ»GKK(I) 

. +AXY»GIJ(I)+AXZ*GIK(I)+AYZ*GJKCI) 

AG = AA(I) / QQR(l) 

Y I = (AG- 1 . ) ♦GSS +AG^DELTA +R(I) 

B = . 5+ ( AQ -1 , ) ♦ ( AXX +AXX +AXY +AXZ) 

Cl = AQ4BXX -Cl. -SIGNX) ♦B 

BI = AQ4BXX -(1. +SIGNX ) ♦» 

AI = -AG ♦ ( BXX +BXX +Q2MBYY +BZZ) ) 

# + ( AG- 1 . ) ♦(Z. ♦(AXX ♦ AYY +AZZ) +AXY +AYZ +AXZ) 

10 RES = ABS ( Y I ) 

IF(RES.LE.FR) GO TO 14 
FR = RES 

IR = I 

JR = J 

KR = K 

14 IF(SIGNX. GT . 0 . ) GO TO 15 

AI = AI +AXT 

Cl = Cl -AXT 

GO TO 16 

15 AI = AI -AXT 

BI = bl + AXT 

16 IF ( SIGN Y .GT . 0 . ) GO TO 17 

AI = AI +AYT 

GO TO 18 

17 AI = AI -AYT 

18 Y I s Y I +AYT^SIGNY^(G(I, JM,K)-GM(I, JM,K)) 

IF(SIGNZ.GT.O. )GO TO 19 

AI = AI +AZT 

GO TO 20 

19 AI s AI -AZT 

20 YI = YI +AZT*S1GNZ* (G(l r J,KM) -GM (I,J,KM)) 

A = l./(AI-BI*C(I-l)) 

C(I) = CI*A 

8 D ( I ) = (YI-8I*D(I-1))»A 

CG =0. 

I =13 

DO 42 M=2,I3 

CG = D ( I ) -C ( I ) *CG 

CORG = ABS ( CG J 

IF (CORG.LE.GD) GO TO 43 
GD = CORG 

IG = I 

JG = J 

KG = K 

43 G(1,J,K) = G(I,J,K) -CG 

42 I =1-1 

J = J -1 

IF ( J-2 ) 61,51,31 

51 IF (ITE2(K) .EG. MX) 13 = LX -1 
GO TO 31 

61 1FCITE2CK) .EQ.MX) GO TO 113 
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IXl = ITE1CK) 

1X2 = 1TE2CK) 

IX1M * 1X1 -1 

IX2P = 1X2 -H 

E = G(1X2,2,K) -G(IX1,2,K) 

DO 100 1= 2 , 1 X 1 M 
M = NX +2 -I 

100 G(I,t,K) = G(M,3,K) -E 

DO 62 1=1X1, 1X2 
XI = A0(1) 

Y 1 = SO ( I , K ) 

X1X1 = XI *X1 

Y1Y1 = Y 1 * Y 1 

HHS = X1X1 +Y1Y1 

DHH = 1. /HHS 

XB = ,5*(X1X1 -Y1Y1) 

XS = XCCK) +XB*SCAL 

X 1 XB = XI *DHH 

X1YB = Y 1 *DtlH 

X1ZBS = -XZCK) * X 1 X B -YZ(K) *X1YB 

Y 1 ZB = XZCK) * X 1 Y B -YZ(K) *XlXb 

GIS = GC1+1,2,K) -G(1-1,2,K) 

GKS s G(l, 2 , K+l ) -G(I,2,K-1) 

GX = AlClJtGIS 

GZ = Cl ( K ) *GKS 

UF = OMEGA*ZS +CAC 

VF = SA 

WF = - (GMEGA*XS +CAS ) 

YXB = -(X1Y8 +X1XB*SX(I,K) ) 

YYb a X 1 XB -XI YB*SX(1,K) 

YZBS = Y 1 ZB -X 1 ZBS*SX ( I , K ) -SZ(I,K) 

X1KSS = X 1 X B * Y X B +XlYb»YYB +X1ZBS*YZBS 

YYSS = YXB+YXB +YYB+YYB +YZBS*YZBS 

HHS = ( U F ♦ Y X B + V F * Y Y B + *IF*YZBS ) *SCAL 

62 G(I,1,K) = GCI,3,K) +CRHS +X1YSS*GX +YZBS*GZ) / ( YYSS*B1 ( 2 ) ) 

DO 102 1= I X 2 P , N X 
M = NX +2 -I 

102 G(i,l,K) = G(M,3,K) + E 
GO TO 103 

113 G ( LX , 2 , K) = G(LX,3,K)*WATY+G(LX-1 ,2,K)*WATX 
DO 114 1= 2, NX 

M = NX+2 -I 

114 GCI,1,K) = G(M,3,K) 

13 = LX -1 

DO 115 1= 2,13 
M = NX+2 -I 

115 G(M,2,K) = G ( 1 , 2 , K ) 

103 CONTINUE 
RETURN 
END 
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APPENDIX D 


LISTING OF SAMPLE DATA 


ROTOR 

BLADES 






FNX 

FNY 

FNZ 





32. 

4, 

6. 





FIT 

covo 

P10 

P20 

P30 

BETAO 

FHALF 

100. 

0.000005 

1.75 

1.0 

1.0 

0.1 

2. 

100, 

0.000005 

1.75 

1.0 

1.0 

0.1 

1 . 

100. 

0.000005 

1.75 

1.0 

1.0 

0.1 

0 . 

FSPEED 

PSI 

ALPHA 

TIPWR 

RAOIUS 

AINF 


80. 

90. 

0.0 

200. 

15. 

334.143 


CREF 

XREF 

FBLADE 

FCLUST 

CDO 



.123 

KPLOTS 

.03075 

1. 

0 . 

0 . 

010101 


FNC 

SWEEP1 

SWEEP? 

SWEEP 

DIHEDl 

DIHED2 

DIHED 

15. 

0.0 

0.0 

0.0 

0 . 

0 . 

0 . 

ZS ( K ) 

XL 

YL 

CHHPD 

THICK 

TWIST 

NEWSFC 

1.00 

-.25 

0.0 

1 . 

1 . 

o.o 

1.0 

YSYM 

FNU 

FNL 





1 . 

65. 

65. 





X 

Y 

UPPER 

SURFACE (M AC A 

00121 




0.000000 

0.000000 

.000250 

.002798 

.000500 

.003945 

.000750 

.004822 

.001000 

.005557 

.001250 

.006203 

.001750 

.007319 

.003000 

.009528 

.003750 

.010622 

.004250 

.011288 

,005750 

.013066 

.006250 

.013603 

.007000 

.014365 

.008250 

.015542 

,009000 

.016202 

.009750 

.016833 

.01 1000 

.017827 

.011750 

.018393 

.012250 

.018759 

.022500 

.024915 

.030000 

.028401 

.037500 

.031374 

.050000 

.035547 

.057500 

.037704 

.062500 

.039027 

.072500 

.041439 

.077500 

.042543 

.082500 

.043587 

.092500 

.045512 

.097500 

.046400 
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.110000 

.048432 

.150000 

.053452 

.170000 

.055287 

.190000 

.056760 

.230000 

.058790 

.250000 

.059412 

.270000 

.059807 

.310000 

.059994 

.330000 

.059820 

.350000 

.059486 

.390000 

.058387 

.410000 

.057643 

.430000 

.056780 

.470000 

.054732 

.490000 

.053560 

.510000 

.052298 

.550000 

.049524 

.570000 

.048021 

.590000 

.046447 

.630000 

.043096 

.650000 

.041325 

.670000 

.039494 

.710000 

.035659 

.730000 

.033658 

.750000 

.031603 

.790000 

.027332 

.810000 

.025117 

.830000 

.022848 

.870000 

.018148 

.890000 

.015715 

.910000 

.013225 

.950000 

.008066 

.970000 

.005393 

.990000 

.002654 

1.000000 

.001260 


ZS(K) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

3.00 

-.25 

0.0 

1 . 

1. 

0.0 

0.0 

ZS ( K ) 

XL. 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

5.00 

-.25 

0.0 

1. 

1 . 

0.0 

0.0 

ZS(K) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

7.00 

-.25 

0.0 

1. 

1 . 

0.0 

0.0 

ZS(K) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

9.00 

- , 2 5 

0.0 

1 . 

1 . 

0.0 

0.0 

ZS ( K ) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

11.0 

-.25 

0.0 

1 . 

1. 

0.0 

0.0 

ZS(K) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

13.0 

- • 25 

0.0 

1. 

1 . 

0.0 

0.0 

ZS(K) 

XL 

YL 

CHORD 

THICK 

TWIST 

NEWSEC 

15.0 

-.25 

0.0 

1. 

1 . 

0.0 

0.0 


%%EOF 
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(a) ONERA straight-tip blade geometry. 



(b) ONERA swept-tip blade geometry. 
Figure 3.- ONERA blade geometry. 



(c) The approximate ONERA straight-tip blade geometry used in the computer code. 



(d) The approximate ONERA swept-tip blade geometry used in the computer code. 

Figure 3.- Concluded. 


52 



ZS = 0.8500 
S Mach = 0.5227 

F 

T 




(b) Advance ratio 0.4 at 30° azimuthal angle. 
Figure 4.- Continued. 


i 
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(c) Advance ratio 0.4 at 60° azimuthal angle. 
Figure 4.- Continued. 
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- 1.2 


ZS = 0.8500 
S Mach = 0.7633 



PSI = 90.0 
F Mach = 0.2406 
T Mach = 0.5976 

ROT17 

• ONERA 




(d) Advance ratio 0.4 at 90° azimuthal angle. 
Figure 4.- Continued. 
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ZS = 0.8500 
S Mach = 0.7311 



PSI = 120.0 
F Mach = 0.2406 
T Mach = 0.5976 

ROT17 

• ONERA 
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(e) Advance ratio 0.4 at 120° azimuthal angle. 
Figure 4.- Continued. 




- 1.2 


ZS = 0.8500 
S Mach = 0.5227 


•8 r 


PSI = 180.0 
F Mach = 0.2406 
T Mach = 0.5976 






PSI = 0.0 
F Mach = 0.3292 
T Mach = 0.5985 


ROT17 

• ONERA 


ZS = 0.9500 
S Mach = 0.5735 






PSI = 30.0 
F Mach = 0.3292 
T Mach = 0.5985 

ROT17 

• ONERA 




i 

| (b) Advance ratio 0.55 at 30° azimuthal angle. 


Figure 5.- Continued. 







ZS = 0.8500 
S Mach = 0.8527 
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(d) Advance ratio 0.55 at 90° azimuthal angle. 
Figure 5.- Continued. 






Figure 5.- Continued. 
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(i) Advance ratio 0.55 at 240° azimuthal angle. 
Figure 5.- Continued. 
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ZS = 0.8500 
S Mach = 0.8192! 


PSI = 60.0 
F Mach =0.3127 
T Mach = 0.6288 
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ZS = 0.9500 
S Mach = 0.8728 



(c) Advance ratio 0.5 at 60° azimuthal angle. 


Figure 6.- Continued. 



ZS = 0.8500 
S Mach = 0.8611 


PSI = 90.0 
F Mach = 0.3127 
T Mach = 0.6288 



Figure 6.- 




Continued . 
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(e) Advance ratio 0.5 at 120° azimuthal angle. 
Figure 6.- Continued. 


76 
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(g) Advance ratio 0.5 at 180° azimuthal angle. 
Figure 6.- Continued. 
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ZS = 0.8500 
S Mach = 0.3920 


PSI = 210.0 
F Mach = 0.3127 
T Mach = 0.6288 




(h) Advance ratio 0.5 at 210° azimuthal angle. 
Figure 6.- Continued. 
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PSI = 270.0 
F Mach = 0.3127 
T Mach = 0.6288 

TFAR1 

• ONERA 



( j ) Advance ratio 0.5 at 270° azimuthal angle. 


Figure 6.- Continued. 





ZS = 0.8500 

S Mach = 0.2775 p SI = 300.0 

F Mach = 0.3127 
9 T Mach = 0.6288 

■ TFAR1 

• • • • ONERA 



(k) Advance ratio 0.5 at 300° azimuthal angle. 


Figure 6.- Continued. 
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